EVSIM  - AN  EVAPORATOR 
SIMULATION  MODEL 
ACCOUNTING  FOR 
REFRIGERANT  AND  ONE 
DIMENSIONAL  AIR 
DISTRIBUTION 


NISTIR  89-4133 


NEW  NIST  PUBLICATION 
August  1989 


ka  ' 

Piotr  A.  Domanski 


U.S.  DEPARTMENT  OF  COMMERCE 
National  Institute  of  Standards 
and  Technology 
Gaithersburg,  MD  20899 


9 


U.S.  DEPARTMENT  OF  COMMERCE 
Robert  A.  Mosbacher,  Secretary 

NATIONAL  INSTITUTE  OF  STANDARDS 
AND  TECHNOLOGY 

Raymond  G.  Kammer,  Acting  Director 


NIST 


NISTIR  89-4133 


EVSIM  - AN  EVAPORATOR 
SIMULATION  MODEL 
ACCOUNTING  FOR 
REFRIGERANT  AND  ONE 
DIMENSIONAL  AIR 
DISTRIBUTION 


Piotr  A.  Domanski 


U.S.  DEPARTMENT  OF  COMMERCE 
National  Institute  of  Standards 
and  Technology 
Gaithersburg,  MD  20899 


Sponsored  by; 

U.S.  Department  of  Energy 
Office  of  Building  and  Community  Systems 
Conservation  and  Renewable  Energy 
Washington,  DC 


August  1989 


U.S.  DEPARTMENT  OF  COMMERCE 
Robert  A.  Mosbacher,  Secretary 

NATIONAL  INSTITUTE  OF  STANDARDS 
AND  TECHNOLOGY 

Raymond  G.  Kammer,  Acting  Director 


ABSTRACT 


This  report  describes  a computer  model,  EVSIM,  of  a refrigerant- to-air  heat 
exchanger  of  the  type  used  in  residential  air  conditioning  as  an  evaporator. 
The  model  provides  performance  predictions  of  a one-slab  or  two-slab 
evaporator  for  a given  refrigerant  enthalpy  at  the  coil  inlet,  saturation 
temperature  and  superheat  at  the  coil  outlet,  and  at  imposed  one  dimensional 
air  mass  flow  distribution  over  the  coil  face. 

The  model  accounts  for  air  distribution  and  for  complex  refrigerant 
circuitry  designs  by  simulating  refrigerant  distribution.  Performance  of 
the  coil  is  calculated  employing  a tube-by-tube  scheme.  Performance  of  each 
tube  is  evaluated  individually  based  on  individual  air  and  refrigerant 
mass  flow  rates  and  their  respective  thermodynamic  states  assigned  for 
each  tube.  The  modelling  effort  emphasis  was  on  the  local  thermodynamic 
phenomena  which  were  described  by  fundamental  heat  transfer  equations  and 
equations  of  state  relationships  among  material  properties. 

This  report  includes  a User's  Guide  and  a listing  written  in  FORTRAN  77. 

Due  to  the  detailed  algorithms  and  tube-by-tube  performance  evaluation  scheme, 
mini  and  main  frame  computers  are  best  suited  for  simulation  studies  using 
EVSIM.  Nevertheless,  the  model  converges  on  an  IBM  AT  compatible  machine 
within  2-6  minutes  when  simulating  a single  slab  evaporator. 


Key  words:  air  conditioner;  coil;  evaporator;  heat  exchanger;  modeling; 
simulation 
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DISCIAIMER 


In  view  of  the  presently  accepted  practice  of  the  building  industry  in  the 
United  States  and  the  structure  of  the  computer  software  used  in  this  project, 
common  U.S.  units  of  measurement  have  been  used  in  this  report.  In  recognition 
of  the  United  States  as  a signatory  to  the  General  Conference  of  Weights  and 
Measures,  which  gave  official  status  to  the  SI  system  of  units  in  1960, 
appropriate  conversion  factors  have  been  provided  in  the  table  below.  The 
reader  interested  in  making  further  use  of  the  coherent  system  of  SI  units 
is  referred  to:  NBS  SP330,  1972  Edition,  'The  International  System  of  Units,' 
or  E380-72,  ASTM  Metric  Practice  Guide  (American  National  Standard  2210.1). 


METRIC  CONVERSION  FACTORS 

Length 

1 inch  (in)  = 25.4  millimeters  (mm) 

1 foot  (ft)  = 0.3048  meter  (m) 

Area 

1 ft^  = 0.092903  m2 

Volume 

1 ft^  = 0.028317  m2 

Temperature 

F = 9/5  C + 32 

Temperature 

Interval 

l^F  = 5/9“C  or  K 

Mass 

1 pound  (lb)  = 0.453592  kilogram  (kg) 

Mass  Per  Unit 
Volume 

1 lb/ft2  = 16.0185  kg/m2 

Energy 

1 Btu  = 1.05506  Kilojoules  (kJ) 

Specific  Heat 

1 Btu/[(lb)(‘’F)]  = 4.1868  kJ/[(kg)(K)] 

Gallon 

1 gallon  = 0.0037854  m^ 
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1. 


INTRODUCTION 


This  report  describes  a computer  model  of  a ref rigerant- to-air  heat  exchanger 
of  the  type  used  as  an  evaporator  in  almost  all  residential  air-conditioning 
applications.  A typical  design  is  shovm  in  Figure  1.  The  refrigerant  flows 
inside  tubes  arranged  in  a certain  circuit  while  the  air  passes  outside  the 
finned  tubes.  The  air  flows  through  the  heat  exchanger,  evaporates  the 
refrigerant,  and  cools  down  in  the  process.  If  the  temperature  of  the  air 
drops  below  the  dew  point  temperature,  water  vapor  contained  in  the  air 
stream  condenses  on  the  heat  exchanger  surface. 


Figure  1.  A schematic  of  a single-slab  evaporator. 
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During  flow  through  the  heat  exchanger  the  temperature  of  the  air  gradually 
decreases  while  the  temperature  of  the  refrigerant  remains  approximately 
constant  as  long  as  the  refrigerant  is  in  the  two-phase  state.  Once  refrigerant 
is  fully  evaporated,  its  temperature  increases  and  the  temperature  difference 
between  the  refrigerant  and  air  decreases.  Air-conditioning  systems  are 
designed  to  operate  at  a small  refrigerant  superheat  at  the  evaporator 
outlet.  Capacity  of  the  evaporator  may  be  significantly  reduced  if  the 
outlet  superheat  is  a result  of  mixing  of  highly  superheated  vapor  and  two- 
phase  flows  leaving  different  circuits  of  the  coil  [1]. 

A number  of  design  features  impact  the  performance  of  an  evaporator  coil. 

Flat,  wavy,  lanced  or  louvered  fins  may  be  used  on  the  air  side.  Tubes  with 
a smooth  inside  surface  are  most  popular  but  enhanced  surfaces  are  also 
available.  Different  materials  may  be  used;  fins  are  normally  made  from 
aluminum  while  tubes  are  manufactured  from  aluminum  or  copper.  Other  * 
significant  design  features  are  staggering  pattern  of  tubes,  tube  diameter, 
number  of  tube  depth  rows,  fin  pitch  and  thickness,  design  of  refrigerant 
circuitry,  and  coil  configuration.  The  last  two  design  aspects  affect 
refrigerant  and  air  distributions. 

Manufacturing  technique  may  also  impact  the  heat  exchanger  performance.  For 
example,  contact  resistance  between  a tube  and  fin  may  depend  on  the  tool 
quality;  outside  surface  treatment  may  promote  better  drainage  of  condensate. 

The  above  description  indicates  that  detailed  modeling  of  an  evaporator 
coil  may  be  quite  involved.  Indeed,  the  complexity  of  all  the  processes 
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associated  with  heat  transfer  between  the  air  and  refrigerant  in  a cross - 
flow,  finned  tube  heat  exchanger  requires  simplifying  assumptions  to  be 
made  during  formulation  of  the  model.  However,  a more  detailed,  fundamentally 
justified  model  has  better  chance  to  correctly  predict  performance  of  the 
modeled  hardware. 

A simple  and  popular  way  to  determine  coil  capacity  depends  on  utilizing 
performance  catalogs  of  major  manufacturers  [2,3],  or  performance  correlations 
developed  by  fitting  the  catalog  data.  Capacity  data  are  presented  in  the 
form  of  charts  for  specific  designs  specified  by  the  tube  pattern,  tube 
diameter,  tube  and  fin  materials  and  the  shape,  spacing  and  thickness  of 
fins.  These  charts  are  usually  developed  for  45®F  refrigerant  saturation 
temperature,  and  80®F  dry  bulb  and  Sy^F  wet  bulb  temperature  of  air.  The 
convenience  and  quickness  of  using  coil  catalogs  comes  at  the  price  of  a 
few  disadvantages;  refrigerant  circuity  can  be  accounted  for  only  in  a 
rudimentary  fashion,  and  extrapolation  from  the  conditions  specified  in  the 
catalog  to  other  operating  conditions  may  be  questionable. 

Computer  based  models  which  provide  performance  predictions  through  evaluation 
of  heat  transfer  relationships  are  more  versatile.  Among  models  in  the 
public  domain,  the  evaporator  model  contained  in  the  general  heat  pump  model 
[4]  was  developed  rather  for  system  studies.  It  uses  effectiveness  vs. 
correlations  and  the  assumption  that  the  heat  exchanger  consists  of  equivalent 
parallel  refrigerant  circuits. 
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Another  evaporator  model,  included  in  the  heat  pump  model  [5],  is  based  on 
a tube-by-tube  approach.  Heat  transfer  in  this  model  is  evaluated  for  each 
tube  independently  based  on  the  refrigerant  temperature  in  the  tube  and  the 
average  air  temperature  for  all  tubes  in  a given  row.  The  selection  of 
tubes  for  performance  evaluation  is  opposite  to  the  refrigerant  flow,  i.e. 
from  the  outlet  to  the  inlet.  This  backward  scheme  results  in  faster 
convergence  of  the  vapor  compression  cycle  simulation  at  the  price  of 
unrealistic  imposition  of  the  same  refrigerant  parameters  at  the  outlet  of 
each  circuit.  The  refrigerant  distribution  in  the  model  is  estimated  based 
on  the  circuitry  layout  and  is  maintained  unchanged  during  coil  simulation. 
The  model  assumes  a uniform  distribution  of  air  and  assigns  the  same  air 
mass  flow  rate  for  each  tube. 

This  report  describes  further  development  of  the  evaporator  model  presented 
in  [5].  The  new  model,  EVSIM,  is  also  based  on  a tube-by-tube  approach  but 
has  several  new  features.  It  can  simulate  performance  of  an  evaporator  coil 
with  non-uniform,  traverse  to  the  tubes,  one -dimensional  air  distribution. 
Refrigerant  distribution  is  simulated  based  on  the  pressure  drops  in  each 
circuit.  A forward  iteration  scheme,  from  the  coil  inlet  to  outlet,  is  used 
which  allows  for  more  realistic  modeling;  refrigerant  superheat  at  different 
circuit  outlets  may  be  different  depending  on  the  heat  gained  and  mass  flow 
in  each  circuit. 
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2. 


MODELING  APPROACH 


2.1  Tube-by-Tube  Method 

The  model  described  in  this  report,  EVSIM,  is  based  on  a tube-by-tube  approach. 
Evaluation  of  performance  for  a single  finned  tube  is  the  basic  part  of  the 
model.  Performance  of  each  tube  is  analyzed  separately  one  at  a time.  Each 
tube  is  associated  with  refrigerant  parameters  and  specific  air  mass  flow 
rate,  inlet  temperature  and  humidity. 

Another  important  part  of  the  model  is  the  logic  which  assigns  tubes  for 
calculation  in  proper  order,  maintains  inventory  of  air  and  refrigerant 
parameters  for  each  tube,  and  maintains  the  convergence  scheme. 

Simulation  starts  with  the  refrigerant  inlet  tube  of  a given  circuit  and 
progresses  consecutively  to  the  following  tubes  until  the  outlet  is  reached. 

If  the  circuit  splits,  the  model  proceeds  with  calculations  of  one  branch 
of  the  circuit  and,  when  the  exit  is  reached,  returns  to  the  split  point  to 
finish  calculations  on  the  remaining  branches.  Once  calculations  for  the 
circuit  are  completed,  the  remaining  circuits  are  calculated  starting  with 
the  inlet  tubes. 

The  selection  scheme  for  tube  calculations  assures  that  refrigerant  parameters 
are  always  known  at  the  entrance  of  the  tube;  they  are  equal  to  the  outlet 
parameters  for  the  proceeding  tube.  At  the  outset  of  simulation  the 
temperature  and  humidity  ratio  of  air  is  known  only  for  tubes  in  the  first 
depth  row  and  is  estimated  for  other  rows.  These  estimates  are  updated  with 
new  calculated  values  as  calculations  progress. 
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2.2 


Air  Distribution 


The  model  can  account  for  non-uniform  air  distribution  between  coil  tubes. 
Each  tube  itself  is  assumed  to  have  uniform  air  distribution  over  its 
length.  Air  distribution  data  have  to  be  provided  to  the  model  in  the  form 
of  values  of  air  velocity  at  the  coil  face  at  discrete  points  on  the 
center  plane,  perpendicular  to  coil  tubes  (specification  of  air  velocity 
at  the  coil  face  is  explained  in  Appendix  A).  From  these  data,  the  model 
derives  velocity  distribution  for  the  face  of  the  coil  and  determines  the 
air  mass  flow  rate  associated  with  each  tube  in  the  first  depth  row. 

Air  mass  flow  rates  associated  with  tubes  past  the  first  row  are  calculated 
based  on  mass  flow  rates  associated  with  the  preceding  tubes.  It  is  assumed 
that  the  general  direction  of  the  air  flow  through  the  slab  is  perpendicular 
to  slab  face  and  that  air  splits  equally  at  each  tube  [1].  Thus  a given 
tube  is  exposed  to  the  air  stream  which  consists  of  50%  of  the  air  streams 
associated  with  the  two  closest  neighbors  in  the  proceeding  row. 

Air  temperature  and  humidity  at  each  tube  past  the  first  row  is  calculated 
by  evaluating  the  mixing  balances  for  air  by  standard  psychometric  equations. 
If  the  mixing  process  results  in  saturated  air,  the  amount  of  condensed 
vapor  is  calculated  and  is  assumed  to  collect  on  fins  and  drain.  The  mixing 
process  is  calculated  by  subroutine  MIXAIR.  The  air  distribution  is  assigned 
by  subroutine  DISTR2. 


6 


2.3  Refrigerant  Distribution 

Refrigerant  distribution  is  governed  in  a coil  by  pressure  drop.  In  this 
report,  separate  consideration  has  been  given  to  single -slab  and  two -slab 
evaporators . 

2.3.1  S ingle - Slab  Evaporator 

EVSIM  can  simulate  one-slab  evaporators  employing  one  expansion  device.  In 
such  a heat  exchanger  refrigerant  pressure  at  each  inlet  tube  is  the  same. 
Refrigerant  mass  flow  rate  through  each  circuit  adjusts  itself  so  pressure 
at  each  outlet  tube  is  also  the  same. 

At  the  outset  of  simulation  the  model  estimates  refrigerant  distribution 
based  on  the  circuitry  layout.  A uniform  resistance  to  flow  at  each  tube  is 
assiimed.  This  assumption  is  crude  since  pressure  drop  will  depend  on  the 
refrigerant  mass  flux  and  heat  transfer  rate.  In  the  course  of  simulation, 
once  pressure  drops  are  calculated  for  each  tube,  the  refrigerant  distribution 
is  updated  using  the  available  pressure  drop  data. 

Figure  2 presents  an  example  of  refrigerant  circuitry  on  one  plane.  Each 
vertical  line  represents  a tube.  Numbers  in  Figure  2 denote  location  of 
each  tube  as  explained  in  Appendix  A. 

The  refrigerant  distribution  in  a coil  is  determined  by  sequential  analysis 
of  each  split  point  and  associated  branches.  The  algorithm  for  evaluation 
of  refrigerant  distribution  was  derived  considering  that  pressure  drop  for 
the  evaporative  flow,  as  evaluated  by  the  Pierre's  correlation  [6],  can  be 
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Refrigerant  Flow 
From  Expansion  Device 

1 


Inlet  Manifold 
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Each  vertical  line  represents  a tube.  A number  next  to  the  line 
location  of  the  tube,  as  explained  in  Appendix  A. 


Figure  2.  An  Example  of  Refrigerant  Circuitry. 
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simplistically  presented  in  the  following  form: 


AP  a f . g2 


(1) 


where  : AP  = pressure  drop 

f = friction  factor 
G = mass  flux 


Considering  that  the  friction  factor  in  the  Pierre's  correlation  is  a 
function  of  the  Reynolds  number  to  the  -0.25  power,  we  can  rearrange 
equation  (1)  to  obtain  the  following  relations  for  each  branch  of  circuitry: 

Ri-Gj-^s^APi  (2) 

where  : R^  = resistance  to  flow  offered  by  a given  branch  leaving  the 

split  point  (it  accounts  for  the  effects  of  tube  geometry, 
fluid  density  and  viscosity) 

= refrigerant  mass  flux  rate  flowing  through  a given 
branch 

AP^  = pressure  drop  in  a given  circuitry  branch 

Equation  (2)  allows  calculation  of  the  flow  resistance  for  each  branch  using 
refrigerant  pressure  drop,  AP^ , and  mass  flux,  , calculated  in  the  previous 
iteration  loop.  Once  flow  resistances,  R^ , are  known,  equation  (2)  allows 
estimation  of  the  ratio  of  mass  fluxes  at  any  two  branches,  considering  the 
fact  that  pressure  drop  through  all  branches  associated  with  a given  split 
point  should  be  the  same. 

Since  the  suun  of  refrigerant  mass  flow  rates  through  the  branches  equals 
the  mass  flow  rate  at  the  split  point,  the  following  equation  holds: 
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1. 


(3) 


n 


where:  = fraction  of  refrigerant  mass  flow  rate  at  the  split  point 

flowing  through  a given  branch 
n = number  of  branches  leaving  the  split  point 

Equation  (3)  and  mass  flux  ratios  generated  by  equation  (2)  constitute  a 
set  of  n equations  with  n unknowns.  This  set  can  be  solved  to  express  the 
fraction  of  refrigerant  mass  flow  rate  at  the  split  point  flowing  through 
any  i*'^  of  its  n branches: 


1. 


(4) 


n 


0.571 


2 (Ri/Rj) 


The  above  equation  is  used  to  update  the  refrigerant  distribution.  As  * 
result  of  distribution  updates,  the  model  iterates  refrigerant  pressure  to 
within  0.05  psi  at  different  circuitry  outlets. 

2.3.2  Two - S lab  Evaporator 

The  model  can  simulate  two  types  of  two-slab  evaporators: 

1.  those  equipped  with  a single  expansion  device  feeding  both  slabs, 

2.  those  equipped  with  two  expansion  devices,  one  for  each  slab. 

These  two  cases  require  somewhat  different  approaches  to  simulate  refrigerant 
distribution.  In  the  first  case,  as  for  a one-slab  coil,  refrigerant 
distribution  is  governed  by  the  refrigerant  pressure  drop  in  the 
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coil  circuits.  Refrigerant  distribution  in  this  case  is  determined  by  EVSIM 
as  for  a one -slab  assembly. 

If  each  coil  of  a two-slab  assembly  is  fed  independently  by  an  individual 
expansion  device,  refrigerant  distribution  is  determined  in  two  stages. 
Firstly,  the  split  of  refrigerant  between  two  slabs  is  determined  solely  by 
the  relative  restrictiveness  of  the  expansion  valves  used.  This  is  due  to 
the  fact  that  expansion  devices  used  in  air-conditioning  applications 
usually  choke  and  exhibit  very  weak  dependence  on  the  evaporator  pressure. 
Secondly,  refrigerant  flow  rate,  already  determined  for  each  slab,  adjusts 
itself  at  individual  circuits  based  on  their  relative  restrictiveness. 

It  should  be  noted  that,  unlike  for  a one  expansion  device  assembly  for 
which  refrigerant  pressures  at  all  inlet  and  outlet  tubes  are  respectively 
equal,  a two  expansion  device  assembly  will  have  the  same  outlet  pressure 
while  the  inlet  pressure  may  be  different  between  the  two  slabs  involved 
(this  will  result  in  different  refrigerant  temperature  profiles  which 
affects  the  heat  transfer  rate) , Since  the  model  is  set  up  to  arrive  at  a 
specified  outlet  pressure  with  the  inlet  pressure  unknown,  an  iterative 
procedure  is  employed  for  the  two  expansion  device  coil. 
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3. 


MODELING  OF  HEAT,  MASS  AND  MOMENTUM  TRANSFER  PROCESSES 


3.1  Heat  Transfer  for  a Tube  in  a Cross-Flow  Arrangement 
If  a single,  separate  tube  of  an  air- to-refrigerant  heat  exchanger  is 
considered,  the  heat  transfer  problem  reduces  to  one  of  a pure  cross-flow. 
For  this  type  of  heat  transfer  closed  expressions  can  be  derived  based  on 
cross-flow  heat  transfer  theory.  According  to  the  general  heat  transfer 
equation: 


Q = U • A . AT„ 


(5) 


where  A = heat  transfer  surface  area 

U = overall  heat  transfer  coefficient 
AT^  = mean  temperature  difference  between  the  heat  exchanging  fluids 
Q = heat  transfer  rate 


In  case  of  a pure  cross-flow  arrangements,  one  of  two  mean  temperature 
differences  applies  for  AT^, , [7]: 


when  temperature  of  one  of  the  (6) 

fluids  does  not  change  (tem- 
perature of  that  fluid  is 
denoted  by  T) 


AT 


m 


^2  ■ ^1 

In  

Ti  - T2  T2 

+ In  — 


^2  ■ 


when  temperature  of  both 
fluids  change 


(7) 


where  T = temperature  of  one  fluid 

t = temperature  of  another  fluid 
AT^j  = mean  temperature  difference 

subscripts  1 and  2 refer  to  tube  inlet  and  outlet 
conditions 
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The  enthalpy  change  of  the  fluid  will  be  according  to  the  equation; 


Q = m • (i2  - ii) 


(8) 


or 


Q = m • Cp  • (T2  - Tj) 


(9) 


Equations  (5)  through  (9)  allow  derivation  of  formulas  for  calculation  of 
the  heat  transfer  rate  in  a tube  with  two-phase  refrigerant  or  superheated 
vapor  flow.  Tubes  in  which  both  flow  patterns  exist  can  be  identified  and 
the  fractions  of  tube  length  associated  with  specific  flow  patterns  can  be 
calculated.  For  refrigerant  pressure  drop  calculations,  the  model  uses 
different  pressure  drop  correlations  for  the  two-phase  portion  and  the 
superheated  vapor  portion  of  the  flow.  Heat  transfer  calculations  employ 
different  inside  tube  heat  transfer  correlations  for  annular  flow  (defined 
by  the  flow  quality  not  greater  than  0.85),  dispersed  flow  (quality  in  the 
range  of  0.85  - 1.00),  and  single-phase,  superheated  vapor  flow  (quality 
greater  than  1.00).  The  derived  correlations  are  presented  below  and  are 
followed  by  explanation  of  the  s3nnbols. 

Annular  Flow 


U • A, 


‘o 


(10) 


Q = ®a  • Op  .(ti  ~ Ti)(l  " exp(~ 


)) 


If  the  calculated  heat  transfer  rate,  Q,  results  in  a refrigerant  quality 
at  the  tube  exit  greater  than  0.85,  the  heat  transfer  rate  within  the  annular 
flow  regime  is  equal  to: 


Q = “r  (ir,85  " ii  ) 


(11) 


13 


The  fraction  of  the  tube  with  flow  quality  up  to  0.85,  ANNUL,  can  be 
calculated  by  the  equation; 


ANNUL 


I»r(ir,8  5 “ , i ) 


U • A„ 


“a  * Cp.aC^i  - Ti)(l  - exp(- 


-)) 


“a  • Cp  3 


(12) 


Dispersed  Flow 


Q = nia 


U 

Cp^^Cl  - ANNUL) (ti  - Ti)(l  - exp( 


(13) 


If  the  calculated  heat  transfer  rate  results  in  the  refrigerant  enthalpy 
exceeding  that  of  the  saturated  vapor,  ij.  y,  the  heat  transfer  rate  within 
the  dispersed  flow  is  equal  to: 


Q = • (ir,v  ~ ir.i) 

The  fraction  of  the  tube  with  dispersed  flow,  XDRY,  can  be  calculated  by 
the  equation: 


* ( Ir  . V , i ) 

XDRY  = (15) 

U . A„ 

I“a  • Cp  .(1  - ANNUL)  (ti  - Ti)(l  - exp( )) 

“a  * Cp,a 

Single-phase  Flow  (superheated  vapor) 

( 1-ANNUL-XDRY ) • Cp  ^ ^ U • A^ 

Q = nij. ‘Op  j.  (t^-T^ ) (1-exp ( ^ (l-exp( )))  (16) 

““r-Cp.r  °*a*Cp,a 
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In  equations  (10)  through  (16)  the  following  nomenclature  was  used: 


Ac 


.P  . r 
1-r  ' i 
, V 
, 85 
“a 
“r 

Ti 

ANNUL 

XDRY 

U 


total  exterior  surface  area  associated  with  the  tube  wetted  by  air 

air  specific  heat  at  constant  pressure 

refrigerant  specific  heat  at  constant  pressure 

refrigerant  enthalpy  at  the  tube  inlet 

enthalpy  of  refrigerant  saturated  vapor 

refrigerant  enthalpy  at  flow  quality  equal  0.85 

air  mass  flow  rate  associated  with  the  tube 

refrigerant  mass  flow  rate  in  the  tube 

air  temperature  upstream  of  the  tube 

refrigerant  temperature  at  tube  inlet 

length  fraction  of  the  tube  with  flow  quality  up  to  0.85 
length  fraction  of  the  tube  with  flow  quality  within  the 
range  0.85  - 1.00 

overall  tube  heat  transfer  coefficient 


3.2  Overall  Heat  Transfer  Coefficient  for  a Finned  Tube 
3.2.1  Dry  Tube 

The  overall  heat  transfer  coefficient,  U,  for  a dry  finned  tube  can  be 
derived  by  summing  up  the  individual  resistances  between  the  refrigerant  and 
the  air , [ 7 ] : 


U = 


r K A„  Xp 

+ + + 

'■  Ap,i  hi  Ap  ,,  kp  \,o 


(Af/Aj(l 


(17) 


where 


Af 

Ao 


Ap  , m 

^ . o 


K 


= fin  surface  area 

= total  exterior  surface  area  exposed  to  air 
= pipe  inside  surface  area 
= pipe  mean  surface  area 
= pipe  outside  surface  area 

= convection  heat  transfer  coefficient  at  the  exterior  surface 
= inside  tube  heat  transfer  coefficient 
= thermal  conductance  of  the  pipe-to-fin  contact 
= thermal  conductivity  of  pipe  material 
= thickness  of  pipe  wall 


T 


f . b 


fin  efficiency 


(18) 


T^  = air  temperature 
Tf  jj  = fin  base  temperature 
Tf  jg  = mean  fin  temperature 
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The  first  term  and  the  fourth  term  of  equation  (17)  refer  to  the  inside  and 
outside  convection  resistances,  respectively.  The  second  term  represents 
the  conductive  heat  transfer  resistance  through  the  tube  wall.  The  third 


term  accounts  for  the  contact  resistance  between  the  outside  tube  surface 
and  the  fin  collar. 

3.2.2  Wet  Tube 

Wet  tube  analysis  is  applicable  to  an  evaporator  when  its  surface  temperature 
is  below  the  dew  point  of  the  air.  As  a result,  moisture  is  removed  from 
the  air  stream  and  by  condensation  on  the  evaporator  external  surface.  At 
evaporator  temperature  above  the  freezing  point  (simulation  range  of  EVSIM) 
the  condensate  flows  down  the  fins  under  the  influence  of  gravity. 

The  heat  transfer  rate  between  the  air  stream  and  the  water  surface  is 
described  by  the  following  equation:  » 


(19) 


The  first  term  accounts  for  sensible  heat  transfer  and  the  second  term 
accounts  for  latent  heat  transfer.  For  air  at  atmospheric  pressure  the 
Lewis  number. 


b, 


(20) 


Le 


is  close  to  1 [7].  Assuming  that  the  fin  efficiency  approximates  the  ratio 
of  the  moisture  content  differences. 


16 


w.  - w. 


(21) 


<t>  = 


f , m 


equation  (19)  assumes  the  following  form  for  a tube  with  flat  fins: 
ifg.w(Wa  - w„)  Af 

dQ  = h,  „(1  + )(1 (1  - .^))(T,  - T^)dA,  (22) 

Cp.aCTa  - T„)  A, 

Symbols  used  in  equations  (19),  (20),  (21),  and  (22)  denote: 

Af  = fin  surface  area 
Aq  = total  external  area 
Cp  a = specific  heat  of  air 

hg  Q = air-side  forced  convection  heat  transfer  coefficient 
hp  p = air-side  mass  transfer  coefficient 
ifg_„  = latent  heat  of  condensation  for  water 
Tg  = temperature  of  air 

= temperature  of  liquid  water  at  the  fin  base 
Q = heat  transfer  rate 
w^  = humidity  ratio  of  air 

= humidity  ratio  of  saturated  air  at  T„  temperature 


The  one -dimensional  heat  conduction  across  the  condensate  film  can  be 
expressed  by  the  equation: 


dQ  = h„.AT,.dA„ 


(23) 


where  h^ 

K 

6 


AT 


W 


K 

• — , heat  transfer  coefficient  for  the  condensate  film 
6 

thermal  conductivity  of  water 
thickness  of  condensate  film 

temperature  difference  across  the  condensate  film 


Using  equations  (19)  and  (23)  and  referring  to  equation  (17),  the  following 
relation  for  the  overall  heat  transfer  coefficient  for  a wet  finned  tube 
can  be  derived: 
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AoXp  1 


U = 


+ + — + 


+ 


^iAp.i  \ ^.o\f 


-1 


1 


(24) 


g , w ('*^a  ) 


)(1 


(l-<^) 


c 


A, 


where  symbols  used  are  defined  as  in  equations  (17)  and  (21). 

The  presented  approach  accounts  for  the  impact  of  moisture  condensation  on 
the  heat  transfer  in  a few  aspects: 

(1)  the  layer  of  condensate  offers  additional  heat  flow  resistance 
(term  3 of  equation  (24)), 

(2)  the  air-side  heat  transfer  resistance  decreases  due  to  effect 

of  condensation  (term  4 of  equation  (24)),  * 

(3)  the  air-side  heat  transfer  coefficient,  h^.  j,,  increases 

since  it  is  sensitive  to  external  surface  geometry  and  the  air  flow 
Reynolds  number  (see  equation  (42)), 

(4)  fin  efficiency  decreases  as  h increases  (see  equation  (36)), 

Equation  (24)  requires  evaluation  of  thickness  of  the  condensate  layer.  The 
analysis  presented  below  allows  evaluation  of  a mean  thickness  of  a condensate 
film  on  a flat,  vertical  fin.  Obviously,  in  a real  application,  condensation 
of  moisture  is  not  uniform  over  the  fin  surface,  and  condensate  thickness 
varies  even  for  flat  fins.  Local  variations  in  the  condensate  thickness  are 
even  more  severe  for  wavy  and  stripped  fins.  At  the  lack  of  a better 
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analytical  approach,  an  average  thickness  of  the  water  layer  is  evaluated 
and  used  in  EVSIM  in  heat  transfer  calculations. 


In  order  to  evaluate  water  layer  thickness,  consider  the  mass  transfer 
equation: 


®«,d  * dw.  = - hu^oCw,  - w^)dAo 


(25) 


For  the  Lewis  number  equal  to  1 equation  (25)  assumes  the  following  form: 


• dw. 


(w„  - w„)  • dA^, 


(26) 


The  change  in  the  air  humidity  ratio  can  be  calculated  by  integrating 
equation  (26),  which  yields: 


Va  . e = . i 


(Wa 


-h. 


- w^)(l  - exp 


p . « 


m 


a , d 


(27) 


Condensation  of  moisture  may  sometimes  occur  only  on  a part  of  the  outside 
surface  associated  with  a tube.  This  may  happen  for  example  for  a tube  in 
which  refrigerant  is  in  a superheated  vapor  state  having  the  inlet  and 
outlet  temperature  below  and  above  the  air  dew  point.  Another  probable  case 
is  when  the  refrigerant  temperature  is  slightly  below  the  air  dew  point.  In 
such  a case  condensation  will  occur  on  the  tube  surface  and  on  that  part  of 
the  finned  area  which  is  below  the  dew  point.  The  fin  surface,  further  from 
the  tube,  having  a temperature  above  the  dew  point  temperature  will  not 
produce  condensation. 
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A rigorous  modeling  of  condensation  on  a part  of  the  outside  surface  requires 
identification  of  surface  areas  above  and  below  the  dew  point.  This  task  is 
extremely  difficult.  The  effect  of  many  variables  affecting  the  temperature 
profiles  of  the  surface  is  unknown.  Notably,  the  effect  of  the  tube 
staggering  pattern  and  temperature  of  the  neighboring  tubes  on  the  fin 
temperature  profile  is  complex.  For  the  above  reasons,  this  part  of  simulation 
is  performed  using  a few  simplifying  assumptions  as  it  was  felt  that  the 
phenomena  are  to  complex  to  model  rigorously  in  a general  evaporator 
simulation  program. 

Moisture  removal  from  the  air  stream  is  calculated  separately  for  a tube 
and  associated  fins.  Fins  associated  with  a given  tube  are  considered  to  be 
circular  and  of  equal  area,  as  shown  in  Figure  3.  A linear  temperature  profile 
of  a tube  surface  between  inlet  and  outlet  is  assumed,  and  portions  of  a 
tube  with  and  without  condensation  are  determined  accordingly.  Also  a ’ 
linear  temperature  profile  is  assumed  for  a fin. 

Mean  temperature  for  a fin  surface,  T^  can  be  expressed  by  the  equation: 

1 


Tf.o  = — jT-dAf 
Af 


(28) 


Applying  a linear  temperature  profile  over  the  fin  and  integrating  we  obtain: 


Tf.B,  = + (Tt  - T„) 

3 2 


(29) 
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Figure  3.  Approximation  method  for  treating  a rectangular-plate  fin 
in  terms  of  a circular-plate  fin  of  equal  area. 


where  = tube  outside  diameter 

D^,  = fin  tip  diameter  (see  Figure  3) 
Tj,  = tube  temperature  at 
Tt  = fin  temperature  at 


Since  the  mean  fin  temperature  can  also  be  related  to  fin  efficiency: 


Tf.m  = T.  - ^(T.  - T,) 


(30) 


the  fin  tip  temperature,  T,.  , and  the  fin  diameter  within  which  condensation 
occurs,  Dg , can  be  determined.  Assuming  that  the  humidity  ratio  of  saturated 
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air  varies  linearly  with  temperature,  the  humidity  ratio  of  the  saturated 
air  corresponding  to  the  mean  temperature  of  the  fin  surface  at  temperature 
below  the  dew  point  can  be  calculated  by  the  following  equation: 


w„  = w„  + (w^ 


(Do 


Do)) 


(31) 


It  should  be  noted  that  EVSIM  approximates  rectangular-plate  fins  by 
circular  plate  fins  just  for  moisture  removal  calculations  only.  The  condensate 
flow  on  fins  and  the  air-side  heat  transfer  calculations  recognize  the  plate 
finned- tube  arrangement. 


The  rate  of  moisture  removal  per  unit  area,  R,  can  now  be  calculated: 


R = °la,d(Va,i  - Va,e)/Ao 


(32) 


where  d ~ mass  flow  rate  of  dry  air 

w^  g = humidity  ratio  of  air  at  tube  row  exit 

Wg  ^ = humidity  ratio  of  air  at  tube  row  inlet 

Assuming  no  air  drag  on  the  liquid  layer,  its  local  velocity  is  expressed 

by  the  closed  solution  of  the  Navier-Stokes  equation  for  a viscous  flow  on 
a vertical  wall: 

Pg 

= — [0.5»y^  - y5]  (33) 

where  = local  liquid  layer  velocity 

p = liquid  density 
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g = gravitational  acceleration 
y = distance  from  the  wall 
8 = liquid  layer  thickness 
H = liquid  absolute  viscosity 


Applying  the  continuity  equation  to  the  liquid  film  of  a unit  width: 


m(2)  = p J dy 

O 


(34) 


and  assuming  uniform  condensation  rate  on  the  fin  (i,e.,  m(z)  = R' • z/h, 
where:  m(z)  = mass  flow  rate  of  condensate  at  elevation  z,  R'  = water 
condensation  rate  by  a fin  of  height  h and  unit  width,  z=  0 at  the  top  and 
z=h  at  the  bottom  of  the  slab) , the  average  condensate  layer  thickness  can 
be  obtained  by  integrating  a local  layer  thickness  over  the  fin  height  and 
dividing  the  obtained  expression  by  the  height.  The  resulting  expression 
is : 


1.082 


.g  * 


(35) 


where  g = gravitational  acceleration 

R'  = condensation  rate  per  unit  width  of  a fin 
= water  dynamic  viscosity 
= water  density 


3.2.3  Fin  Efficiency 

The  addition  of  fins  to  the  tubes  greatly  increases  the  outer  heat  transfer 
area  but  at  the  expense  of  decreasing  the  mean  temperature  difference 
between  the  surface  and  the  air  stream.  The  parameter  called  fin  efficiency, 
4>,  is  used  to  rate  the  thermal  effectiveness  of  a fin. 
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Fin  efficiency  for  a circular  flat  fin  on  a singular  tube  was  discussed  by 
Gardner  in  his  classic  paper  [8].  Gardner  solved  the  differential  equation 
for  describing  the  temperature  distribution  and  presented  fin  efficiency 
curves . 


In  the  case  of  residential  air-source  evaporators,  fins  are  not  circular  but 
continuous,  rectangular  plates  serving  all  tubes  in  the  slab.  The  shape  of 
fin  serving  a particular  tube  depends  on  the  tube  pitch  and  transverse  tube 
spacing  in  the  assembly.  The  method  proposed  by  Schmidt  [9]  and  described  in 
[10]  is  sensitive  to  the  pattern  of  tube  staggering  and  is  employed  here. 

The  fin  efficiency,  4> , is  calculated  in  terms  of  the  fin  root  radius,  r^  , 
and  two  parameters,  M and  6: 


M 


'2  . h ] 
. • t , 


(36) 


e = (R/r„  - 1)(1  + 0.35-ln(R/r„)  (37) 

tanh  (M«rp»tf) 

4 = (38) 


where  h = 


r ifg.w(Wa  - w^)  1 

.o  1 + 

'■  Cp,.  (T,  - T,)  J 


if  w^  > w^ 


, o 


otherwise 


R = equivalent  fin  tip  radius 
r^  = outside  radius  of  tube 
t = fin  thickness 
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The  staggered  tube  configuration  results  in  a hexangular  fin  as  shown  in 
Figure  4.  The  two  dimensions  shown  in  the  figure  can  be  defined  interchan- 
geably. They  must  be  defined  in  such  a way  to  have  Y always  greater  than 


Figure  4.  Definition  of  dimensions  for  calculating  fin  efficiency 
by  Schmidt's  method. 

or  equal  to  y.  Following  the  notation  shown  in  the  figure,  the  ratio  of  the 
equivalent  fin  tip  radius  and  the  fin  root  radius  is  calculated  by  the 
following  equations: 


^ = y/^o 

^ = Y/y 

R/r  = 1.27*V(/3  - 0.3)°-5 


(39) 

(40) 

(41) 
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Numerous  assumptions  are  associated  with  the  calculation  of  the  fin 
efficiency.  The  fin  is  asstimed  to  be  thin  and  the  air-side  heat  transfer 
coefficient  to  be  constant  over  the  fin  surface.  Several  practical 
aspects  of  the  heat  exchanger  operation  are  not  taken  into  account  due  to 
inability  to  model  and/or  perception  that  their  impact  on  the  fin  efficiency 
value  is  insignificant.  For  example,  heat  transfer  between  neighboring 
fins  is  not  taken  into  account.  Similarly,  the  effect  of  the  discontinuities 
in  the  fin  material  (lanced  fins)  on  the  fin  efficiency  is  ignored  in  the 
model;  it  will  vary  from  one  lance  design  to  another  and  is  generally  unknown. 
More  study  is  needed  to  include  the  above  effects. 

3.3  Forced  Convection  Heat  Transfer  at  the  Air-Side  of  a Dry  Plate 
Finned-Tube 

EVSIM  is  concerned  with  modeling  of  evaporators  having  continuous  fins  on 
a staggered  array  of  circular  tubes.  Three  consecutive  sections  describe  air- 
side  heat  transfer  correlations  applicable  to  three  different  fin  designs. 
Calculation  of  the  air-side  heat  transfer  coefficient  is  handled  by  subroutine 
AIRHT3 . 

3.3.1  Flat  Fins 

The  correlation  of  Gray  and  Webb  [11]  was  selected  to  calculate  the  air- 
side  heat  transfer  coefficient  for  flat  fins.  This  correlation  was  developed 
using  laboratory  data  on  16  heat  exchangers  applying  a multiple  regression 
technique.  The  rms  error  was  7.3%. 

The  correlation  provides  an  average  value  for  the  j -factor  for  a heat  exchanger 
with  four  or  more  tube  depth  rows  (no  change  in  the  j -factor  after  4-rows 
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is  assumed).  It  has  the  following  form: 


= 0.14*Re 


-0.32 


fs  1 

-0.502 

f \ 

s 

Si. 

Po. 

0.0312 


(42) 


where  j, 


hj.  o • Pr^  / 3 


Gp  ' C„  ^ 


j -factor  for  four  or  greater  number  (43) 

of  depth  rows 


hg  o = outside  surface  forced  convection  heat  transfer  coefficient 
Pr  = Prandtl  number 

Gg  = air  mass  flux  based  in  the  minimum  flow  area 
Cp  g = specific  heat  of  air  at  constant  pressure 
D,  . G, 

Re  = , Re3molds  number 

Dp  = outside  diameter  of  tube 

H = dynamic  viscosity 

= tube  spacing  normal  to  air  flow 

= tube  spacing  in  air  flow  direction  (depth  row  pitch) 
s = spacing  between  adjacent  fins 


To  calculate  an  average  value  for  the  j -factor  for  heat  exchangers  with  less 
than  four  depth  rows,  jj,  (where  N<4)  , Gray  and  Webb  [11]  provided  the  following 
equation: 


U =U  • 0.991[2.24-Re"°-°®2(N/4)-o.03i jo.607(4-N) 


(44) 


where  = value  obtained  by  equation  (42) 

Tube-by-tube  simulation  requires  the  availability  of  the  air-side  heat 
transfer  coefficient  for  a tube  in  a given  row.  Assuming  that  each  row 
weights  equally  on  the  average  air-side  heat  transfer  coefficient  of  the 
coil,  the  heat  transfer  coefficient  value  for  the  depth  row  N,  jj,  ^ , can 
be  approximated  by  the  formula: 
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Jn,R  “ 


(45) 


where  , Ju-i  = average  j -factors  for  heat  exchangers  with  N and 

N-1  depth  rows,  respectively,  obtained  by  equation 
(43)  or  (44) 


3.3.2  Wavy  Fins 

Comprehensive  data  on  performance  of  flat  and  wavy  fins  were  published  by 
Beecher  and  Fagan  [12].  Webb  and  Trauger  [13]  used  their  results  as  a base 
for  their  correlation  which  are  applied  in  EVSIM  to  provide  the  value  of 
the  enhancement  of  the  air-side  heat  transfer  coefficient  for  a wavy  fin 
over  a flat  fin.  This  enhancement  value  is  then  used  as  a multiplier  to  the 
heat  transfer  coefficient  calculated  for  a flat  fin  by  equations  presented 
in  Section  3.3.1 

Beecher  and  Fagan  reported  fins  performance  in  terms  of  the  Nusselt  number 
based  on  the  arithmetic  mean  temperature  difference  between  the  air  and 
refrigerant,  Nu^j.  . Correlations  provided  by  Webb  and  Trauger  express  this 
Nusselt  number  with  the  Graetz  number,  Gz,  and  non-dimensionalized  geometric 
parameters.  For  wavy  fins  the  correlations  are: 


for  Gz  < 25 


Nu,r  = 0.5*Gz°®® 


fs  1 

0.11 

f 

s 

-0.09 

rs  1 

0.12 

r2Spi 

Do- 

Dc. 

Si. 

.Si  . 

-0.34 


(46) 


for  Gz  > 25 


Nu,r  = 0.83-Gz°-^® 


fs  1 

0.13 

f 

S 

-0.16 

fs  1 

0.25 

f2Spl 

Dc. 

Dc. 

Si. 

Si  . 

-0.43 


(47) 
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These  equations  predict  88%  of  the  base  data  within  ±5%,  and  99%  of  the 
data  within  10%. 

Correlations  for  flat  fins  have  the  following  form: 


for  Gz  < 25 


f % 

s 


Nu^j.  = 0.4*Gz°  • 


vDc 


-0.23 


N 


0.23 


(48) 


for  Gz  > 25 


= 0.53-Gz°®2 


-0.23 


N 


0.31 


(49) 


Flat  fin  correlations  predict  98%  of  the  source  data  [12]  within  5% 


To  convert  Nu^j.  to  the  Nusselt  number  based  on  the  logarithmic  mean  temperature 
difference,  Nu,  the  following  equation  is  used  [12]: 


1 + 2-Nu.r/Gz 

Nu  = 0.25»Gz*ln  (50) 

1 - 2.Nu„/Gz 


The  symbols  used  in  equations  (46)  through  (50)  are  defined  as  follows: 


h . Dh 

Nu  = , Nusselt  number  (51) 

k 

Gz  = Re»Pr  , Graetz  number  (52) 

N*Si 
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2*s(l  - P) 


, hydraulic  diameter 


(53) 


- 

7-(l  - )3)  + 2-s.^/D, 
P*V,.Dh 

Re  = , Rejmolds  number 

M 

»r*D^ 

^ = 

4.S,.Si 

7 = [1  + 4(Sd/2.Sp)2]0-5 
Vc=  V„„/(l  - 


(54) 


(55) 


(56) 

(57) 


where  N = number  of  tube  depth  rows 

^max  ~ velocity  based  on  minimum  flow  area 

Sjj  = fin  pattern  depth,  tip  to  valley  (see  Figure  5) 
Sp  = fin  pattern  length,  tip  to  tip 
s = spacing  between  adjacent  fins 
p = density 

k = thermal  conductivity 
H = d3mamic  viscosity 


The  variables  characterizing  the  wave  pattern  were  set  within  subroutine 
AIRHT3  as  follows:  Sp  = 0.5  • , 

Sj  = s 

Dg  = Dp  + 2*t  where  t = fin  thickness. 


3.3.3  Lanced  Fins 

Lanced  fins  are  those  enhanced  fins  which  have  arrays  of  small  strips 
raised  from  the  base  plate.  Nakayama  and  Xu  [14]  proposed  a heat  transfer 
correlation  for  such  fins.  Their  formula  is  in  the  form  of  a heat  transfer 
correlation  for  a flat  fin  and  a multiplier  which  provides  correction  for 
heat  transfer  enhancement  due  to  the  raised  strips. 
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Figure  5,  Geometry  of  a wavy  fin. 
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To  obtain  a common  simulation  base  for  different  heat  transfer  surfaces, 
EVSIM  uses  this  multiplier  in  conjunction  with  the  equations  presented  for 
a flat  fin  in  Section  3.3.1. 


The  lanced  fin  enhancement  multiplier  is  a function  of  the  geometry  parameters 
shown  in  Figure  6.  The  proposed  correlation  has  the  following  form  [14]: 


"t 


Fj  = 1 + 1093 


1 . 2 A 


fi-  2 . 0 9 


+ 1.097 


2 . 2 6 T}  . 6 6 


(58) 


strip  area  (2nj. -l)!^ 

where  = = (59) 

total  fin  area  St»Sj-  0.25»?r»Do^ 

was  set  to  0.275  in  subroutine  AIRHT3.  The  applicable 
range  of  <f>^  is  0.2  - 0.35. 

P*Vn>ax-De 

Re  = , Reynolds  number  (60) 

A* 

Dg  = hydraulic  diameter  of  the  minimum  free  flow  area 

Ug  = number  of  strips  in  the  enhanced  zone 

Ig  = width  of  a strip 

Sg  = length  of  a strip 

t = fin  thickness 

s = spacing  between  adjacent  fins 

^max  ~ velocity  in  the  minimum  flow  area 
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Figure  6.  Geometry  of  a lanced  fin. 


3.4  Forced  Convection  Heat  Transfer  Inside  an  Externally 
Heated  Tube 

Tubes  with  smooth  inner  surface  are  predominantly  used  in  residential  heat 
exchangers,  however,  tubes  with  enhanced  surfaces  are  also  available.  A 
significant  amount  of  research  has  been  performed  over  the  years  on  heat 
transfer  for  smooth  surfaces.  Since  a smooth  surface  is  clearly  defined, 
laboratory  data  of  different  researchers  can  be  pooled  together  for  the 
correlation  development.  As  a result,  a number  of  correlations  for  the  heat 
transfer  coefficient  inside  a smooth  tube  are  available. 
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The  situation  is  different  for  tubes  with  enhanced  surfaces.  Much  less 
research  have  been  done  on  enhanced  tubes  than  on  smooth  tubes,  and  there 
is  an  infinite  number  of  ways  the  inner  surface  may  be  enhanced.  There  is 
not  adequate  information  in  the  open  literature  to  allow  calculation  of  the 
enhanced  tube  heat  transfer  coefficient  in  a generalized  way.  Consequently, 
while  the  internal  tube  heat  transfer  coefficient  for  smooth  tubes  is 
calculated  in  EVSIM  by  a well  founded  correlation,  the  heat  transfer 
coefficient  for  enhanced  tubes  is  evaluated  by  multiplying  the  heat 
transfer  values  for  smooth  tubes  by  the  fixed  correction  factors.  The 
following  discussion  on  flow  patterns  refers  to  smooth  tubes  and  may  not 
necessarily  be  accurate  in  all  aspects  for  tubes  with  enhanced  surfaces. 

Refrigerant  enters  the  evaporator  from  an  expansion  device  at  a quality  of 
about  20%  and  forms  annular  flow  almost  instantly.  The  quality  increases  as 
the  flow  proceeds  and  the  liquid  layer  at  the  wall  gets  thinner.  At  higher 
qualities,  reported  by  different  researchers  to  be  in  the  range  from  65%  to 
95%,  the  refrigerant  vapor  has  enough  kinetic  energy  to  gradually  destroy 
the  liquid  layer.  Consequently,  the  annular  flow  is  followed  by  the  mist 
flow  and  single-phase  superheated  vapor  flow. 

3.4.1  Single-Phase  Flow 

Smooth  tubes . The  single -phase  forced  convection  heat  transfer  coefficient, 
hgp,  for  a smooth,  heated  tube  is  calculated  by  a commonly  used  equation 
[15]: 

h^p  = 0.023*Re°®-Pr°-^-k/Di  (60) 
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where  Re  = Reynolds  number 
Pr  = Prandtl  number 

k = thermal  conductivity  of  refrigerant  vapor 

= tube  inner  diameter 

Enhanced  tubes . The  heat  transfer  coefficient  for  enhanced  inner  surfaces 
is  calculated  within  EVSIM  by  multiplying  h^^^  obtained  from  equations  (60) 
by  a correction  factor  equal  to  2.0.  The  value  of  the  selected  correction 
factor  is  an  average  of  the  heat  transfer  enhancement  reported  by  Khanpara 
et  al.  [16]  for  micro-fin  tubes  with  R-22. 

3.4.2  Two -Phase  Flow  with  Evaporation 

Refrigerant  flow  with  evaporation  is  subdivided  in  the  model  in  two  flow 
patterns;  annular  flow  and  mist  flow.  The  quality  value  of  0.85  was  selected 
in  the  model  as  the  border  point  between  these  two  flow  patterns. 

Smooth  tubes,  annular  flow.  A correlation  developed  by  Gungor  and  Winterton 

[17]  is  used  in  EVSIM  to  calculate  the  evaporative  heat  transfer  coefficient 
for  the  annular  flow  regime  in  smooth  tubes.  This  correlation  was  developed 
with  the  aid  of  a large  data  bank  which  included  4300  points  by  28  authors 
covering  7 fluids.  The  form  of  the  correlation  is  consistent  with  Chen's 

[18]  approach  in  that  it  recognizes  two  distinct  mechanisms  for  the  heat 
transfer;  nucleate  boiling  and  forced  convection.  The  two-phase  evaporation 
heat  transfer  coefficient,  h^^ , is  expressed  as  a weighted  average  of  the 
convective  single-phase  heat  transfer  coefficient,  h^ , and  the  pool  boiling 
heat  transfer  coefficient,  hpo^^ , responsible  in  the  correlation  for  the 
nucleate  boiling  contribution  to  the  heat  transfer: 


(61) 


\n  = E’hi  + S*hpooi 
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hi  = 0.023*Rei°  • ® ‘Pr/ • ^ ‘ki/Di  (62) 

\ooi  = 55-P/-12(_iog^^p^)-0.55.M-0.5.q0.67  (w/m^ -K)  (63) 

E = 1 + 24000-Boll®  + 1.37.X-0-86  (64) 

S = (1  + 1.15-10~®-E2-Reii^)-i  (65) 

In  the  case  of  a horizontal  tube  and  the  Froude  number,  Fr,  smaller  than 
0.05,  E and  S should  be  multiplied  by  E2  and  $2,  respectively: 

E2  = Fr(° • i-2-Fr ) (66) 

S2  = Fr°-®  (67) 

$ 

The  symbols  used  in  equations  (61)  through  (67)  have  the  following  meaning: 

G(1  - x)Di 

Re^  = , liquid  Reynolds  number  (68) 

Ml 

G2 

Fr  = , Froude  number  (69) 

-g 

Di  = tube  inside  diameter 
G = refrigerant  mass  flux 
g = acceleration  of  gravity 
M = molecular  weight 
Py  = reduced  pressure 
Pr^  = liquid  Prandtl  number 
q = heat  flux  (W/m^ ) 

X = flow  quality 

= liquid  thermal  conductivity 
Pj  = liquid  density 
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Bo 


Q 


, boiling  number 


Smooth  tubes,  mist  flow.  The  heat  transfer  coefficient  for  the  mist  flow, 
hn,  (flow  quality  range  0.85  - 1.00),  is  calculated  in  EVSIM  by  weighting  by 
quality  the  heat  transfer  coefficient  values  obtained  by  equation  (60)  and 
equation  (61): 

h„  =(1.0-x)h,„  + (x-0.85)h,p  (70) 


where  x = average  fractional  flow  quality  for  the  mist  flow  in  a tube 
hgj,  = heat  transfer  coefficient  obtained  by  equation  (61)  for 
flow  quality  0.85 

hgp  = heat  transfer  coefficient  obtained  by  equation  (60)  for 
saturated  vapor  flow 


Enhanced  tubes.  The  heat  transfer  coefficient  for  an  enhanced  inner 
surface  is  calculated  by  applying  a multiplier  of  1.45  to  the  heat  transfer 
coefficient  calculated  for  a smooth  tube  within  the  two -phase  region.  The 
selected  multiplier  value  is  an  average  of  the  enhancement  range  (1.2  to 
1.7)  reported  by  Khanpara  et  al.  [16]  for  a micro-fin  tube  with  R-22. 


3.5  Thermal  Conductance  of  the  Pipe-to-Fin  Contact 

The  thermal  conductance  of  the  pipe-to-fin  contact  is  calculated  by  the 
correlation  provided  by  Sheffield  et.  al.  [19]. 


htf  = exp 


'' 

rl  • FPI  • d-|0.75 

1.251 

6.902  + 2.889  • 

..  1 

t • FPI 

(71) 
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where : 


= thermal  conductance  (Btu/h»ft^  •?) 

I = tube  expansion  interference  (in) 

FPI  = number  fins  per  inch  (1/in) 

d = indentation  diameter  of  Vickers  microhardness  test, 
25g  load  (/im)  (d  is  set  in  EVSIM  to  44.2  fim) 

Dp  = tube  outside  diameter  (in) 

t = fin  thickness  (in) 


This  correlation  is  restricted  to  mechanically  expanded  copper  tubes  with 
aluminum  fins.  Other  applicability  limits  which  has  to  be  satisfied  are: 


Dp  < 0.625  (in) 

FPI  < 18 

0.003  < I < 0.0079  (in) 

3.6  Refrigerant  Pressure  Drop  in  a Tube 

The  total  pressure  drop  for  flow  in  a tube  consists  of  the  pressure  drops 
due  to  friction,  momentum  change  and  gravity.  The  gravitational  pressure  drop 
is  neglected  in  EVSIM. 

Comments  made  in  the  beginning  of  chapter  3.4  about  availability  of  heat 
transfer  performance  data  for  enhanced  surfaces  are  applicable  also  for 
pressure  drop.  Consequently,  pressure  drop  for  enhanced  surfaces  is  calculated 
applying  correlations  for  smooth  tubes  and  predetermined  correction  factors. 

3.6.1  S ingle - Phase  Flow 

Smooth  tubes.  Frictional  pressure  drop  can  be  calculated  by  the  Fanning 
equation  with  the  Fanning  friction  factor  as  per  equations  (72)  and  (73): 
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dP  2•f^G2 


(72) 


dL  ‘p 


f = 0.046*Re‘° -2 


(73) 


Pressure  drop  due  to  momentum  change  can  be  calculated  by  the  following 
equation: 


dP  dv 

— = - G2  — (74) 

dL  dL 


where  P = pressure 

L = coordinate  along  the  tube  axis 
G = refrigerant  mass  flux 
= tube  inner  diameter 
V = refrigerant  specific  volume 


Enhanced  tubes.  Frictional  pressure  drop  is  calculated  by  equations  (72) 
and  (73)  applying  multiplier  of  1.5  to  the  obtained  result.  This  multiplier 
provides  a half  of  the  enhancement  provided  in  EVSIM  for  the  respective 
inside  tube  heat  transfer  coefficient. 

3.6.2  Two-Phase  Pressure  Drop  with  Evaporation 

Smooth  tubes . The  pressure  drop  for  two-phase  flow  with  evaporation  can  be 
calculated  by  the  correlation  proposed  by  Pierre  [6] , based  on  his  experiments 
with  refrigerants  R-12  and  R-22.  Pierre's  correlation  combines  the  frictional 
and  momentum  change  effects  in  one  equation: 
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(75) 


L Ax 

AP  = (f  — + — )G2.v„ 

X 


where  f = friction  factor  calculated  by  equation  (76) 
Xjj,  = mean  quality 
Ax  = quality  change 

= Vl  + (Vy  “ Vl  ) , mean  specific  volume 


f = 0.0185 


.25 


iReJ 


(76) 


J*if6*Ax*gc 

where  = (77) 

L*g 

G*Di 

Re  = (78) 

/^L 

J = 778.17  (Ibf ‘ft/Btu) , mechanical  equivalent  of  heat 
g =32.2  (ft/s^),  gravitational  acceleration 

gg  = 32.2  (ft* lb/(lbf. ) , dimensional  constant 


The  formula  for  the  friction  factor  contains  the  Reynolds  number  and  the 
term,  Kj , referred  to  by  Pierre  as  a boiling  number,  making  the  friction 
factors  sensitive  to  vapor  generation  rate  at  the  liquid-vapor  interface. 
Pierre's  correlation  was  verified  in  [20]  providing  overall  better  agreement 
with  experimental  data  for  pressure  drop  of  refrigerants  R-12  and  R-22  than 
the  other  most  popular  correlation,  that  of  Martinelli  and  Nelson  [21]. 
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Enhanced  tubes . Pressure  drop  for  tubes  with  enhanced  surfaced  is  calculated 
applying  a multiplication  factor  1.225  to  the  pressure  drop  value  obtained 
from  equation  (75).  The  factor  1.225  represents  a half  of  enhancement 
provided  in  EVSIM  for  the  two-phase  heat  transfer  coefficient  for  these 
surfaces . 
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4. 


MODEL  VERIFICATION 


The  verification  of  the  model  consisted  of  comparing  the  coil  capacity 
prediction  with  test  data  at  different  air  distribution.  Laboratory  test 
data  of  an  evaporator  coil  at  four  different  air  velocity  profiles  presented 
in  [1]  were  used  for  model  verification.  The  velocity  profiles  projected  on 
the  coil  face  were  used  as  input  to  the  model  along  with  other  operating 
conditions  reported  for  these  tests.  Test  conditions  and  laboratory  and 
simulation  results  are  shown  in  Table  1.  The  air  velocity  profiles  are 
presented  in  Figure  7. 

Figure  7 illustrates  that  the  measured  air  velocity  profiles  differ  from  a 
parabolic  profile  associated  with  fully  developed  flow  in  a straight  duct. 

The  author  believes  that  the  odd  profiles  are  due  to  the  duct  configuration, 
short  entrance  length  and  the  presence  of  the  coil  in  the  duct.  Various  air 
velocity  profiles  were  obtained  by  changing  the  coil  angle  with  respect  to 
the  duct.  It  should  be  noted  that  the  velocity  profiles  at  different  coil 
angles  were  the  result  not  of  the  coil  position  alone  but  rather  of  the 
combined  effects  of  change  of  the  distance  between  coil  surface  and  the  duct 
entrance  and  the  duct  bend  past  the  coil. 

Total  capacity  results  are  also  shown  in  Figure  8.  The  change  of  capacity 
is  presented  as  a function  of  the  coil  angle.  The  figure  displays  the  model 
sensitivity  to  the  air  distribution.  Predictions  of  total  capacities  are 
within  8.2  percent  of  the  test  results.  The  model  predicted  degradation  trend 
for  65°  and  45°  angles  except  at  25°  where  the  prediction  of  capacity  is 
within  1 percent  of  the  test  result.  A number  of  reasons  could  cause  the 
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Table  1.  Results  of  Laboratory  Tests  and  Simulation  Runs. 
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VELOCITY  (tt/s)  VELOCITY  (ft/s) 


a = 90"  a = 65" 


a = coil  angle  as  shown  in  Figure  8 
location  = distance  from  slab  edge  shown  as  X coordinate  in  Figure  A3 


Figure  7.  Air  velocity  profiles. 
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CAPACITY  (kBtu/h) 


a (deg) 


Figure  8.  Simulation  and  test  results. 
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break  in  the  prediction  trend  between  25°  and  other  angles.  One  possible 
reason  may  be  the  fact  that,  due  to  lack  of  space  needed  for  the  Pitot  tube 
used  for  velocity  measurement,  the  air  velocity  profiles  for  the  first  three 
angles  were  obtained  from  measurements  past  the  coil  while  the  25°  velocity 
profile  was  measured  before  the  coil  face.  It  is  also  possible  that  the 
velocity  profiles  derived  from  the  air  velocity  measurements  at  the  center 
plane  represented  at  different  degree  the  average  air  flow  over  the  length 
of  each  tube  at  different  coil  angles. 

Prediction  of  the  latent  capacity  is  not  as  good  as  of  the  total  capacity. 
The  highest  discrepancy  between  test  and  simulation  results  is  19.2%.  Lack 
of  a closer  agreement  is  not  a complete  surprise  since  the  mass  transfer  and 
the  condensate  flow  are  very  difficult  to  simulate  and  are  represented 
simplistically  in  EVSIM.  This  simulation  aspect  requires  further  study  to 
allow  for  simulation  improvement.  ’ 
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5. 


CONCLUDING  REMARKS 


The  evaporator  model  presented  in  this  report,  EVSIM,  possesses  two  useful 
features : 

- ability  to  account  for  air  distribution 

- ability  to  simulate  refrigerant  distribution  in  complicated  circuits 
These  two  capabilities  are  in  a way  compatible  since  refrigerant  distribution 
in  a coil  depends  not  only  on  circuitry  specification;  change  of  air 
distribution  affects  heat  transfer  impacting  refrigerant  quality  and 
pressure  drop  thus  changing  distribution  of  refrigerant.  As  a result  of 

it's  complexity,  the  model  is  more  likely  to  provide  satisfactory  capacity 
predictions  than  would  simpler  models,  particulary  in  cases  where  air 
distribution  is  not  uniform  over  a coil  face  or  refrigerant  circuitry  is 
complex.  The  model  should  be  very  helpful  as  a coil  design  tool, 
decreasing/eliminating  the  need  for  development  tests,  since  detailed 
information  on  performance  of  each  tube  and  refrigerant  superheat  at  each 
circuit  outlet  can  be  obtained  from  the  model. 

Practical  aspects  of  using  the  model  include  preparation  of  a coil  data 
file  and  computer  requirements.  The  detailed  performance  information  which 
EVSIM  can  provide  comes  at  the  price  of  the  need  for  a detailed  input  data 
file  of  the  heat  exchanger.  In  addition  to  general  design  input  like  tube 
diameter,  fin  spacing  and  thickness,  tube  staggering  pattern,  CFM,  etc., 
specification  of  refrigerant  circuitry  and  air  distribution  is  needed. 

The  model  is  written  in  FORTRAN  77.  Due  to  detailed  algorithms  and  tube-by- 
tube  performance  evaluation  scheme,  mini  and  main  frame  computers  are  best 
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suited  for  simulation  studies  using  EVSIM;  however  the  model  has  been 
installed  on  an  IBM  AT  compatible  computer  converging  within  2-6  minutes  for 
a single  slab  coil,  and  up  to  14  minutes  for  a two  slab  coil. 

The  model  is  capable  of  simulating  the  performance  of  evaporators  with 
flat,  wavy  or  lanced  fins  on  the  air  side,  and  tubes  having  smooth  or 
enhanced  inner  surfaces.  Capacity  predictions  for  flat  fin  and  smooth  tube 
coils  are  expected  to  be  most  reliable  since  these  basic,  plain  surfaces  are 
well  defined  and  their  respective  heat  transfer  correlations  are  based  on 
extensive  data  bases.  Capacity  prediction  for  evaporators  with  enhanced 
surfaces  may  not  be  as  accurate  since  each  enhanced  surface  is  different 
with  its  own  heat  transfer  characteristics;  general  correlations  may  not 
well  represent  their  performance.  This  is  particulary  true  for  inner  tube 
surfaces  for  which  no  general  correlations  are  available  and  simple  enhancement 
multipliers  are  used  to  account  for  performance  change  with  respect  to  a * 
smooth  tube.  Other  area  in  which  new  algorithms  could  be  applied,  if 
available,  to  enhance  the  model  are  condensation  of  water  vapor  and  prediction 
of  the  coil  latent  capacity. 
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APPENDIX  A.  PROGRAM  USER'S  GUIDE 


The  simulation  program  of  an  air-source  evaporator,  EVSIM,  predicts  performance 
of  a given  evaporator  at  imposed  air  conditions  and  CFM,  and  refrigerant 
parameters  at  the  heat  exchanger  inlet  and  outlet.  The  model  is  based  on  a 
tube -by- tube  approach  and  account  for  non-uniform  distributions  of  air  and 
refrigerant . 

The  program  is  written  in  Fortran  77  and  makes  use  of  standard  Fortran 
mathematical  functions.  It  consists  of  a main  module,  evaporator  simulation 
subroutine,  EVPHX2 , and  30  subprograms.  Capability  of  the  model  and 
relations  used  for  fluid  mechanics,  heat  transfer  and  mass  transfer 
calculations  are  described  in  the  main  text  of  this  report.  Equations  used 
in  moist  air,  water,  and  refrigerant  property  routines  are  described  in  [5]. 

Al . Input  Data 

Input  data  are  read  by  the  program  from  data  files  and  a terminal  (batch 
file)  depending  on  data  category. 

Al . 1 Data  Read  from  Data  Files. 

To  run  the  program,  a user  must  establish  two  data  files  on  the  system, 

DATAREF  and  DTEV,  for  refrigerant  constants  and  evaporator  data. 

Refrigerant  property  constants.  Constants  for  evaluation  of  thermodynamic 
and  transport  properties  for  Refrigerant  12  and  Refrigerant  22  are  shown  in 
Table  Al  and  Table  A2 , respectively.  Constants  of  the  selected  refrigerant 
have  to  be  in  a disk  file  named  DATAREF  from  which  they  are  read  by  subroutine 
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DATAIN.  Comments  statements  inserted  into  subroutine  DATAIN  help  to 
identify  specific  constants. 

Evaporator  Data.  Evaporator  data  are  read  by  subroutine  RDATA3  from  a file 
named  DTEV.  Coding  of  evaporator  data  in  DTEV  is  described  in  Table  A3 
which  includes  Fortran  symbols  and  their  short  explanation.  The  convention 
by  which  DTEV  has  to  be  prepared  and  some  of  the  variable  names  are 
explained  in  detail  below. 

EVSIM  can  simulate  one-slab  and  two-slab  evaporator  assemblies  (coils).  The 
input  format  shown  in  Table  A3  is  complete  for  a one-slab  coil;  a data  file 
for  a two-slab  assembly  will  have  lines  3 through  25  repeated  as  lines  26 
through  48  with  input  data  for  slab  # 2.  An  example  of  a data  file  for  the 
two-slab  assembly  of  Figure  Al  is  shown  in  Table  A4. 

For  two-slab  coils,  the  model  can  handle  coils  in  which  each  slab  is  fed  by 
an  individual  expansion  device,  or  both  slabs  are  fed  by  a single,  common 
expansion  device.  These  arrangements  are  coded  in  the  data  file  in  line  2. 
In  the  present  version,  EVSIM  is  unable  to  handle  a slab  fed  by  more  than 
one  expansion  device. 

The  last  input  datum  at  line  7,  SFLOW,  is  the  expansion  device  scaling 
factor.  For  capillary  tubes,  it  may  be  determined  with  the  aid  of  Fig.  39 
of  [22].  This  scaling  factor  was  included  into  a data  file  to  allow  EVSIM 
to  calculate  the  refrigerant  split  between  two  slabs  in  two-slab  assemblies 
if  separate  expansion  devices  are  used  for  each  slab.  For  evaporator 
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assemblies  employing  just  one  expansion  device  (both  single  and  two-slab 
coils)  the  value  of  SFLOW  has  no  bearing  on  calculations;  any  real  number 
may  be  inputed  to  satisfy  the  format  of  the  READ  statement.  For  evaporators 
with  two  expansion  devices,  the  ratio  of  SFLOW(l)  and  SFL0W(2)  has  to  be 
equal  to  the  ratio  of  refrigerant  mass  flow  rates  through  expansion 
devices  associated  with  respective  slabs.  The  ratio  of  flow  factors,  as 
presented  in  Fig.  39  of  [22],  satisfies  this  condition;  hence  flow  factors 
are  prescribed  as  input  to  the  evaporator  data  file.  In  the  most  probable 
case,  the  two  expansion  devices  are  identical;  in  such  situation  two 
identical  real  numbers,  say  1.,  satisfy  the  format  of  the  READ  statement 
and  modeling  requirements  of  EVSIM. 

Line  3 contains  general  dimensions  of  a slab.  The  s3mibols  used  are  explained 
graphically  in  Figure  A2.  Note  that  BSPACE  denotes  the  distance  between  the 
edge  of  the  open  channel  for  air  flow  to  the  center  of  tube  #1 . Numbering 
of  tubes  should  be  done  as  shown  in  Figure  Al  and  Figure  A2  for  two -slab 
and  one-slab  coils,  respectively.  Numbering  should  start  at  the  leftmost  tube 
in  the  first  tube  row  (facing  the  incoming  air).  Once  the  last  tube  of  a 
given  row  is  given  a number,  a consecutive  number  is  given  to  the  leftmost 
tube  in  the  next  depth  row.  It  is  not  important  at  which  coil  side  tube 
numbers  are  assigned  so  long  as  it  is  consistent  with  the  later  assignment 
of  the  air  velocity  profile  for  the  coil. 

Line  5 contains  tube  information.  ISUR  is  an  identifier  for  smooth  or 
enhanced  inner  surface.  Only  for  smooth  surface  are  the  heat  transfer  and 
pressure  drop  calculations  performed  using  well  establish  correlations. 
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Because  of  the  variety  of  enhanced  surface  and  lack  of  confirmed  correlations, 
the  heat  transfer  and  pressure  drop  for  enhanced  surfaces  are  evaluated  by 
applying  correction  factors  to  smooth  tube  calculations  (refer  to  section 
3 . 4 and  3.6). 

Line  6 groups  air  side  fin  data.  Heat  transfer  for  plain,  wavy  and  lanced 
fins  is  calculated  by  separate,  dedicated  correlations.  Wave  pattern  for  wavy 
fins  and  strip  set  for  lanced  fins  are  assigned  within  subroutine  AIRHT3 
(see  section  3.3). 

Lines  8 through  20  describe  refrigerant  circuitry.  Description  of  refrigerant 
circuitry  depends  on  specification  for  every  tube  the  tube  which  supplies 
it  with  refrigerant.  It  is  done  in  the  numerical  order  starting  from  the  tube 
numbered  as  1.  in  line  8,  with  ten  tubes  per  line.  Taking  as  an  example  the 
circuitry  shown  in  Figure  Al , since  the  first  field  of  line  8 is  designated 
for  tube  #1,  tube  # 2 shall  be  placed  in  this  field  since  tube  #2  feeds 
tube  #1.  If  a given  tube  receives  refrigerant  from  an  inlet  manifold,  the 
input  shall  be  zero.  EVSIM  can  handle  coils  with  up  to  130  tubes  per  slab. 
Enter  999  in  the  data  field  if  the  tube  does  not  exist. 

Line  22  through  25  are  for  air  velocity  measurement  data  input.  The  data 
consist  of  location  of  the  measurement  and  the  air  velocity  at  this  location. 
The  measurements  should  be  taken  in  the  central  plane,  perpendicular  to 
coil  tubes.  Based  on  these  measurements,  EVSIM  develops  the  air  velocity 
profile  in  the  central  plane  and  uses  it  to  evaluate  the  air  mass  flow 
rate  for  individual  tubes  in  the  assembly.  Figure  A3  shows  how  EVSIM 
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utilizes  the  air  velocity  data.  A velocity  profile  is  created  by  straight 
line  interpolation  between  the  discrete  data  points.  Note  that  the  air 
velocity  at  the  sides  of  the  slab  is  assigned  equal  to  velocities  of 
respective  leftmost  and  rightmost  data  points  (VX(2,1)  and  VX(2,6)  in  Figure 
A3).  The  maximum  data  input  is  16  points.  The  minimum  data  input  is  one 
point.  The  latter  case  results  in  prescription  of  a uniform  velocity 
profile  over  the  slab  face. 

It  should  be  noted  that  air  velocity  data  serve  to  establish  the  velocity 
profile  only  and  are  not  used  as  input  of  air  volumetric  flow  through  the 
coil.  The  velocity  profile  is  integrated  to  obtain  a CFM  value  which  is  then 
compared  by  the  program  with  CFMTOT  to  prorate  local  velocities  so  CFM  and 
CFMTOT  are  equal. 


Al . 2 Data  Read  from  a Terminal 


Data  read  from  a terminal  (batch  file)  Include: 

- air  parameters : dry  bulb  temperature 

relative  humidity 

- refrigerant  parameters:  inlet  quality 

outlet  saturation  temperature 
outlet  superheat 

- program  output  controlling  parameters. 


This  input  is  solicited  by  the  program  during  program  execution  in  the 


following  format: 


1.  Request: 
Response : 


DATE: 

alphanumeric  response  up  to  16  characters 


2.  Request: 


IPR  = 0 FOR  MAIN  RESULTS  OUTPUT  ON  THE  DEFAULT  DEVICE 
IPR  = 1 FOR  MAIN  RESULTS  OUTPUT  TO  FILE  "RESULT" 

IDIA  = 0 FOR  NO  ADDITIONAL,  DETAILED  RESULTS  OUTPUT 
IDIA  = 1 FOR  ADDITIONAL,  DETAILED  RESULTS  OUTPUT 
TOGETHER  WITH  MAIN  RESULTS 

IDIA  = 2 FOR  ADDITIONAL,  DETAILED  RESULTS  OUTPUT  TO  FILES 
"DIAG"  AND  "DIAGO" 
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IPR,  IDIA  = 

Response:  two  integer  numbers  separated  by  a comma 

3.  Request:  TAIR  = AIR  DRY  BULB  TEMPERATURE,  (F) 

RH  = AIR  RELATIVE  HUMIDITY,  (DECIMAL  FRACTION) 

TAIR,  RH  = 

Response : two  real  numbers  separated  by  a comma 

4.  Request:  XI  = REFRIGERANT  INLET  QUALITY,  (DECIMAL  FRACTION) 

TSAT2=  REFRIGERANT  SAT.  TEMP.  AT  COIL  OUTLET,  (F) 

TSUP2=  REFRIGERANT  SUPERHEAT  AT  THE  COIL  OUTLET,  (F) 

XI,  TSAT2,  TSUP2  = 

Response : three  real  numbers  separated  by  commas 

A2 . Output  Data 

EVSIM  concludes  a run  when  it  converges  on  the  imposed  refrigerant  saturation 
temperature  and  superheat  at  the  evaporator  outlet  with  the  following 
convergence  parameters: 

saturation  temperature:  ± 0.05  ®F 
superheat:  ± 2.0  ®F. 

Once  the  model  converges  within  2.0  °F  of  superheat,  the  intermediate  results 
from  the  last  two  iteration  loops  are  used  to  interpolate  the  performance 
results  to  the  superheat  value  specified  in  the  input  data.  If  the  model 
is  unable  to  converge,  it  still  performs  interpolation  and  provides  a 
warning  message . 

A short  results  output  and  a detailed  results  output  are  available  from 
EVSIM,  as  indicated  in  the  previous  section.  The  short  results  version 
contains  a short  summary  of  input  data  and  coil  performance  information  from 
the  last  iteration  loop  which  include  state  information  on  refrigerant 
leaving  the  individual  outlet  tubes.  The  last  part  of  the  output  is  the  coil 
performance  at  the  requested  superheat  at  the  coil  exit.  An  example  of  the 
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short  results  output  is  given  in  Appendix  E. 

A more  detailed  results  output  (IDIA=2)  contains  the  intermediate  results 
for  individual  iteration  loops  (file  DIAG)  and  counters  indicating  the 
number  of  iteration  loops  required  to  converge  in  a given  run  (IDIAO) . The 
option  of  a detailed  output  may  be  used  if  more  specific  information  about 
the  coil  performance  is  sought;  this  output  may  be  redesigned  to  suit  an 
individual  need. 
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Table  Al . Property  Constants  for  Refrigerant  12  in  the  Input  Format  to  Program 
EVSIM 


REFRIGERANT  12 

6 . 9330000E+02 , 5 . 9690000E+02 , 2 . 8700000E-02 
4 . 5967000E+02 . 1 . 8505300E-01 , 2 . 7182818E+00 
9 . 1835883E+01 , -7 . 9131381E+03 , -1 . 2471522E+01 
1.0892245E-02,0. ,0. 

120. ,312. 

3 . 4840000E+01 , 1 . 8691368E+01 , 2 . 1983963E+01 
5 . 3341175E+01 , -3 . 1509939E+00 , 1 . 

5 . OOOOOOOE-01 , 3 . 3333333E-01 . 2 . 

8 . 8734000E-02 , 6 . 5093890E-03 . 0 . 

-3 .4097270E+00 , 1 . 5943480E-03 , -5 . 6762767E+01 
6.0239450E-02. -1.8796180E-05,1.3113990E+00 
-5.4873700E-04,0. ,0. 

0. ,3.4688340E-09,-2.5439070E-05 
0.,0..0. 

0. ,5.4750000E+00 

8 . 0945000E-03 , 3 . 3266200E-04 , -2 . 4138960E-07 
6.7236300E-11,0. ,0. 

3 . 9556551E+01 , -1 . 6537940E-02 

0.75800014, -0.44230204E-02,0.22659166E-04 

-0.80936502E-07,0.48639626E-09,-0.38992915E-ll 

0.16672029E+01,-0.37690766E-01,0.48997637E-03 

-0.32598077E-05,0.10701433E-07,-0.14068541E-10 

0.0262, 5. 8E-05,0. 

0. ,0. ,0. 

-0.10065153E01,0.37970707E-01,-0.53504633E-03 

0.36446277E-05,-0.12031578E-07,0.15478212E-10 

0.49000002E-01,-0.11950213E-03,0.36320252E-07 

0.52080296E-09,-0.31689972E-ll,-0.31688932E-13 

0. 33482605E00, -0. 10634881E-01 ,0. 14902322E-03 

-0.10213785E-05,0.34027736E-08,-0.44299642E-ll 

0.0043,1.7E-05,0. 

0. ,0. ,0. 

-0.12814194E00,0.48691398E-02,-0.68452384E-04 

0.46711199E-06,-0.15483721E-08,0.20037774E-ll 

0.21699998E00,0.14159610E-03,0.64704511E-06 

0.55390004E-08,-0.13778530E-10,-0.17912061E-12 

-0. 15409625E01 ,0.64801723E-01 , -0 . 91255037E-03 

0.62140834E-05,-0.20422284E-07,0.26255064E-10 
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Table  A2 . Property  Constants  For  Refrigerant  22  in  the  Input  Format  to  Program 
EVSIM 


REFRIGERANT  22 

6 . 6450000E+02 , 7 . 2191000E+02 , 3 . 0525000E-02 
4.5967000E+02,1.8505300E-01,2.7182818E+00 
6 . 7598246E+01 , -8 . 8538843E+03 , -7 . 8610310E+00 
5 . 0448235E-03 , 4 . 4574700E-01 , 6 . 8610000E+02 
140. ,240. 

3 . 2760000E+01 , 5 . 4634409E+01 , 3 . 6748920E+01 
-2.2292560E+01,2.0473289E+01,3.3333333E-01 
6.6666667E-01 ,1. ,1.3333333E+00 
1 . 2409800E-01 , 2 . OOOOOOOE-03 , 0 . 
-4.3535470E+00,2.4072520E-03,-4.4066868E+01 
-1 . 7464000E-02 , 7 . 6278900E-05 , 1 .4837630E+00 
2.3101420E-03,-3.6057230E-06,0. 

-3 . 7240440E-05 , 5 . 3554650E-08 , -1 . 8450510E-04 
1 . 3633870E+08 , - 1 . 6726120E+05 , 0 . 

5 . 4820000E+02 , 4 . 2000000E+00 

2.8128360E-02,2.2554080E-04,-6.5096070E-08 

0. ,0. ,2.5734100E+02 

6 . 2400900E+01 , -4 . 5333500E-02 

0.64600E+00,-0.29194E-02,0.12164E-04 

-0.74985E-07,0.83951E-09,-0.37512E-ll 

0.69684E01.-0.24319E00,0.35924E-02 

-0.26187E-04,0.93884E-07,-0.13301E-09 

0.26600E-01,0.63804E-04,0.10761E-06 

-0.32061E-08,0.43463E-10.-0.13175E-12 

-0.66330E00,0.25757E-01,-0.37913E-03 

0.27734E-05,-0.10053E-07,0.14474E-10 

0.63000E-01,-0.15820E-03,-0.12289E-06 

0. 17453E-08, -0.52685E-11 . -0.10857E-13 

-0.46705E00,0.19421E-01,-0.28507E-03 

0.20429E-05,-0.71963E-08,0.99460E-ll 

0.48000E-02,0.19881E-04,0.24815E-08 

0.28518E-09,-0.62001E-ll,0.31001E-13 

0.51539E00,-0.18601E-01,0.26762E-03 

-0.18936E-05,0.6589lE-08,-0.90041E-ll 

0 . 27100E+00 , 0 . 24054E-03 , 0 . 38936E-07 

0.23481E-07,-0.97345E-10,0.44953E-12 

0.49002E00,-0.83123E-02,0.13105E-03 

-0.96884E-06,0,36462E-08,-0.52089E-ll 
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Table  A3.  Input  Data  Code  for  an  Evaporator  Assembly 


The  input  data  described  below  constitute  a complete  data  set  for  a single 
slab  evaporator.  All  input  data  are  in  Fortran  free  field  input  format 
with  data  values  on  the  same  line  separated  by  commas. 


Line  1 : COILID 

COILID  = alphanumeric  coil  information,  maximum  70  characters 


Line  2 : NSLABS , 
NSLABS 

I EXP 

CFMTOT 


I EXP,  CFMTOT 

= number  of  heat  exchanger  slabs  in  the  coil  assembly, 
possible  values:  1 or  2 

= number  of  expansion  devices  in  the  assembly,  possible 
values:  1 if  NSLABS  = 1 

1 or  2 if  NSLABS  = 2 

= total  volumetric  flow  of  air  through  the  assembly, 

(f t^ /min) 


Line  3 : SLABID 

SLABID  = alphanumeric  slab  information,  maximum  70  characters 


Line  4 : 


BSIDE(l),  BSPACE(l),  WIDTH(l)  (see  Figure  A2) 

BSIDE(l)  = height  of  the  coil,  (in)  * 

BSPACE(l)  = distance  between  the  edge  of  the  coil  and  location  of 
tube  # 1,  (in) 

WIDTH(l)  = width  of  a coil,  equal  to  the  length  of  tubes  exposed 
to  the  duct  air,  (in) 


Line  5 : TPCH(l) , 
TPCH(l) 
DPCH(l) 

Did) 

D0(1) 

TKM(l) 

ISUR(l) 


DPCH(l),  DI(1),  D0(1),  TMK(l),  ISUR(l)  (see  Figure  A2) 

= tube  pitch  in  each  depth  row,  (in) 

= distance  between  neighboring  tube  depth  rows,  see  Fig. 
Al,  (in) 

= tube  inside  diameter,  for  grooved  tubes  use  the  minimum 
diameter,  (in) 

= tube  outside  diameter,  (in) 

= thermal  conductivity  of  a tube  material,  (Btu/(ft 'h* F) ) 
= 1 for  a smooth  inner  surface 
= 2 for  an  enhanced  inner  surface 


Line  6 : FPCH(l), 
FPCH(l) 
FTK(l) 
FMK(l) 
IFIN(l) 


FTK(l),  FMK(l),  IFIN(l) 

= center  to  center  distance  between  fins,  (in) 

= fin  thickness,  (in) 

= fin  material  thermal  conductivity,  ( (Btu/ft*h*F) ) 
= 1 for  flat  fins 
= 2 for  wavy  fins 
= 3 for  lanced  fins 
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Line  7 : NTUB(1,1) , 

NTUB(1 , 1) 

NTUB(1,2) 
NTUB(1 ,3) 
NTUBd  ,4) 
NTUB(1,5) 
SFLOW(l) 


NTUB(1,2),  NTUB(1,3).  NTUB(1,4),  NTUB(1,5),  SFLOW(l) 
= number  of  tubes  in  the  first  depth  row  (facing  the 
incoming  air) 

= number  of  tubes  in  the  second  depth  row 

= number  of  tubes  in  the  third  depth  row 

= number  of  tubes  in  the  forth  depth  row 

= number  of  tubes  in  the  fifth  depth  row 

= expansion  device  flow  factor 


Line  8 : IFR0M(1,I),  I = 1.10 

IFR0M(1,1)  = number  of  the  tube  from  which  tube  1 receives  refrig. 
IFR0M(1,2)  = number  of  the  tube  from  which  tube  2 receives  refrig. 
IFR0M(1,3) 


IFR0M(1,9) 

IFROM(1,10) 

Line  9 : IFR0M(1,I). 

IFROM(l,ll) 

IFR0M(1,12) 

IFR0M(1,13) 


number  of  the  tube 
I = 11,20 

number  of  the  tube 
number  of  the  tube 


from  which  tube  10 


from  which  tube  11 
from  which  tube  12 


receives  refrig 

receives  refrig 
receives  refrig 


Line 

10: 

Line 

11: 

Line 

12: 

Line 

13: 

Line 

14: 

Line 

15: 

Line 

16: 

Line 

17: 

Line 

18: 

Line 

19: 

Line 

20: 

Line 

21: 

IFR0M(1,19) 

IFROM(1,20) 

IFR0M(1,I) , 

IFR0M(1,I) , 

IFR0M(1,I) , 

IFR0M(1,I) , 

IFR0M(1,I) , 

IFR0M(1,I) . 

IFR0M(1,I) , 

IFR0M(1,I) , 

IFR0M(1,I) , 

IFR0M(1,I) , 

IFR0M(1,I) , 

NTEST(l) 

NTEST(l) 


= number  of  the  tube  from  which  tube  20  receives 
I = 21,30 
I = 31,40 
I = 41,50 
I = 51,60 
I = 61,70 
I = 71,80 
I = 81,90 
I = 91,100 
I = 101,110 
I = 111,120 
I = 121,130 


= number  of  air  velocity  measurement  points, 
possible  values:  minimum  2,  maximum  16 


refrig 
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Line  22:  Xfl.N')  . N=1 , 8 

X(l,l)  = location  of  the  first  air  velocity  measurement  point 

(distance  between  the  edge  of  the  slab  closest 
to  tube  #1  and  the  velocity  measuring  probe,  see 
Figure  A3),  (in) 

X(l,2)  = location  of  the  second  air  velocity  measurement  point 

X(l,3)  = location  of  the  third  air  velocity  measurement  point, 

if  non-existent,  input  0.0,  (in) 


X(l,8)  = location  of  the  eight  air  velocity  measurement  point, 

if  non-existent,  input  0.0,  (in) 


Line  22:  X(1,N) , 
X(l,9) 

X(l,10) 


N=9,16 

= location  of  the  ninth  air  velocity  measurement  point, 
if  non-existent,  input  0.0,  (in) 


X(l,16)  = location  of  the  sixteenth  air  velocity  measurement 

point,  if  non-existent,  input  0.0,  (in) 


Line  23:  VX(l.N) . N=l,8 

VX(1,1)  = air  velocity  at  the  first  measurement  point,  (ft/s) 

VX(1,2)  = air  velocity  at  the  second  measurement  point,  (ft/s) 

VX(1,3)  = air  velocity  at  the  third  measurement  point, 

if  non-existent,  input  0.0,  (ft/s) 

VX(1,8)  = air  velocity  at  the  eight  measurement  point, 

if  non-existent,  input  0.0,  (ft/s) 


Line  25:  VX(N) . N=9,16 

VX(1,9)  = air  velocity  at  the  ninth  measurement  point,  (ft/s) 

VX(1,10) 


VX(1,16)  = air  velocity  at  the  eight  measurement  point,  (ft/s) 

if  non-existent,  input  0.0,  (ft/s) 


Line  25  completes  data  file  for  a one  slab  evaporator  (slant  coil).  If  two 
slab  evaporator  is  considered  (V-shape,  A-shape  coil),  the  data  file  has  to 
contain  additional  23  lines  of  data  in  which  the  second  slab  is  described. 
In  this  case,  lines  3 through  48  are  dedicated  to  slab  # 1,  and  identical 
lines  26  through  48  describe  slab  # 2.  Subroutine  RDATA3 , which  reads  the 
evaporator  data  file,  assigns  index  '2'  instead  of  '1'  when  reading  lines 
26  through  48  (e.g.,  reads  VX(2,16)  instead  VX(1,16)).  Table  A4  provides  an 
example  of  a data  file  for  an  A-shape  coil. 
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Table  A4.  Example  of  a Data  File  for  an  A-shape  Coil 
Below  is  the  coding  of  the  coil  shown  in  Figure  Al . 


***DTEV***  a- shape  coil,  3 DEPTH  ROWS,  16  TUBES  PER  ROW. 

2,2,1120.0 

***DATA  FOR  SLAB  # 1*** 

16.0625.0. 8125.17.875 

1.00. 0.875.0.363.0.394.223. .1 
0.0789,0.008,128. ,2 

16.16.16.0. 0.1. 

2,3,19,5,6,22,23,7,8,25 

10,27,12,30,14,15,33,17,18,4 

37.21.39.0. 9.11.28.45.13.31 
32,48,34,35,36,20,38,41,40,24 

42.25.26.43.44.29.46.47.999.999 

999.999.999.999.999.999.999.999.999.999 

999,999,999,999,999,999,999,999,999,999 

999,999,999,999,999,999,999,999,999,999 

999,999,999,999,999,999,999,999,999,999 

999,999,999,999,999,999,999,999,999,999 

999,999,999,999,999,999,999,999,999,999 

999,999,999,999,999,999,999,999,999,999 

999.999.999.999.999.999.999.999.999.999 
6 

1.0.  4.. 7.. 10.. 12. 2. 14. 8. 0.0. 0.0 

0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. 

3. 3. 4. 8. 6. 1.4. 5. 4. 2. 4. 4.0.  .0. 

0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. 

***DATA  FOR  SLAB  # 2*** 

16.0625.0. 8125.17.875 

1.00. 0.875.0.363.0.394.223. .1 
0.0789,0.008,128. ,2 

16.16.16.0. 0.1. 

2,3,19,5,6,22,23,7,8,25 

10,27,12,30,14,15,33,17,18,4 

37.21.39.0. 9.11.28.45.13.31 
32,48,34,35,36,20,38,41,40,24 

42.25.26.43.44.29.46.47. 999 . 999 

999.999.999  999  999  999,999,999,999,999 

999.999.999.999.999.999.999.999.999.999 

999,999,999,999,999,999,999,999,999,999 

999,999,999,999,999,999,999,999,999,999 

999,999,999,999,999,999,999,999,999,999 

999,999,999,999,999,999,999,999,999,999 

999,999,999,999,999,999,999,999,999,999 

999,999,999,999,999,999,999,999,999,999 
6 

1.0.  4.. 7.. 10.. 12. 2. 14. 8. 0.0. 0.0 
0.,0.,0.,0.,0.,0.,0.,0. 

4. 1.4. 9. 4. 3. 2. 9. 3. 9. 3. 3.0.  .0. 

0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. 
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Tubing  bends  on  visible  side 


Figure  Al . Example  of  circuitry  specification  for  an  A-shape  coil. 
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Tubing  bends  on  visible  side 
Tubing  bends  on  the  other  side 


Figure  A2.  Specification  of  slab  circuitry  and  dimensions. 
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Air 

Flow 


Figure  A3.  Example  of  air  velocity  measurement  points  and  velocity  profile. 


APPENDIX  B.  Listing  of  the  Program,  EVSIM 

The  list  of  functions  and  subroutines  used  in  EVSIM  is  given  in  Table  Bl. 
Figures  Bl  and  B2  present  flow  charts  of  the  main  program  and  the  evaporator 
simulation  subroutine,  respectively.  The  following  pages  contain  listing  of 
the  main  program,  and  subroutines  and  functions  in  the  alphabetical  order. 
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Table  B1 . Functions  and  Subroutines  used  in  EVSIM 


NAME 


AIRPR 

AIRHT3 

BALFLl 

CPCV 


DATAIN 

DISTR2 

DYNADP 

EVGUNG 

EVPDP2 

EVPHX2 

FINCON 

FINEF2 

FRACT 

HYDDIA 

ITRPR2 

MIXAIR 

OVLWET 

PHWET2 

RDATA3 

SATP 

SATPR 

SATT 

SATVF 

SPHDP 

SPHDPl 

SPHTC 

TRACE3 

VPSV 

VPVHS 

WATPR 


PURPOSE 


Calculate  properties  of  wet  air 

Calculate  air  side  heat  transfer  coefficient  for  flat,  wavy 
or  lanced  fin 

Adjust  refrigerant  flow  distribution  in  an  evaporator  based 
on  pressure  drop  in  different  circuits 

Calculate  specific  heat  at  constant  volume  and  at  constant 
pressure,  specific  heat  ratio,  and  sonic  velocity  of  refrigerant 
vapor  from  temperature  and  pressure 
Read  file  DATAREF  with  refrigerant  constants 
Determine  air  distribution  for  each  tube  of  the  coil 
Calculate  dynamic  pressure  drop  for  a flow  in  a tube 

Calculate  forced-convection  evaporative  heat  transfer  coefficient 
Calculate  frictional  evaporation  pressure  drop  in  a tube 
Simulate  performance  of  an  evaporator  coil 

Calculate  the  thermal  conductance  for  a fin- to -fin  contact 
Calculate  fin  efficiency 

Calculate  refrigerant  distribution  in  a split  point  based 
on  pressure  drop  in  the  downstream  circuits 
Calculate  hydraulic  diameter  for  air  flow  through  a slab 
Calculate  refrigerant  vapor  thermodynamic  properties  from  given 
pressure  and  specific  volume  or  enthalpy  or  entropy 
Calculate  properties  of  the  air  stream  resulted  from 
the  mixing  process  of  two  wet  air  streams  > 

Calculate  overall  heat  transfer  coefficient  for  a wet 
finned  tube 

Calculate  refrigerant  thermodynamic  properties  from  given  pressure 
and  enthalpy 

Read  file  DTEV  with  evaporator  coil  data  and  prepare 
these  data  for  EVPHX2 

Calculate  refrigerant  saturation  pressure  from  given  temperature 
Calculate  refrigerant  dynamic  viscosity  and  thermal  conductivity 
for  liquid  and  vapor,  and  liquid  specific  heat  at  saturation 

Calculate  refrigerant  saturation  temperature  from  given  pressure 
Calculate  specific  volume  of  saturated  liquid  refrigerant  from 
given  temperature 

Calculate  pressure  drop  for  a single-phase  flow  in  a smooth  tube 
Calculate  pressure  drop  for  a single-phase  flow  in  a smooth  tube 
Calculate  force-convection  single-phase  heat  transfer  coefficient 
in  a smooth  tube 

Estimate  refrigerant  distribution  in  an  evaporator  coil  based  on 
circuitry  configuration 

Calculate  specific  volume  of  vapor  from  given  temperature  and 
pressure 

Calculate  refrigerant  vapor  thermodynamic  properties  from  given 

temperature  and  pressure 

Calculate  water  and  frost  properties 
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start 


j 


INPUT: 

Coil  Data 

Refrigerant  Constants 
Air  Conditions 
X1,TSATMSUP" 


I 

- 

Estimate  PI 

,X1.RMS 

1 

— 

Simulate  Evaporator 

Symbols 


TSAT* 


Yes 

Print  Results 

Stop 


inlet  pressure 
refrig,  mass  flow  rate 
calculated  saturation 
temperature  at  coil 
outlet 

requested  saturation 
temperature  at  coil 
outlet 

calculated  superheat 
at  coil  outlet 
requested  superheat 
at  coil  outlet 
inlet  quality 


Figure  Bl . Logic  of  the  main  program  of  EVSIM. 
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^ San  j 


INPUT: 

Col  Geometry 
kiM  Air  Condition 
Ratiig.  Outlet  Condition 
Retiig.  Inlet  Quality 


Symbols 


- refrig,  enthalpy  obtained 
from  present  iteration 
loop 

- refrig,  enthalpy  obtained 
from  previous  iteration 
loop 

- refrig,  enthalpy  at 
coil  exit  obtained  from 
present  iteration  loop 

- refrig,  enthalpy  at 
coil  exit  obtained  from 
previous  iteration  loop 


Figure  B2 . Logic  of  the  evaporator  simulation  subroutine,  EVPHX2 . 
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PROGRAM  EVSIM 


••••  THIS  IS  THE  MAIN  PROGRAM  OF  THE  EVAPORATOR  MODEL  — EVSIM  — 

THIS  PROGRAM  IS  A PART  OF  THE  REPORT: 

"EVSIM  - AN  EVAPORATOR  SIMULATION  MODEL  ACCOUNTING  FOR  REFRIGERANT 
AND  ONE  DIMENSIONAL  AIR  DISTRIBUTION" 

PREPARED  AT  NATIONAL  INSTITUTE  OF  STANDARDS  AND  TECHNOLOGY 
AUTHOR:  P.A.  DOMANSKI 

NIST  PRINCIPAL  INVESTIGATOR:  D.A.  DIDION 
SPONSOR:  U.S  DEPARTMENT  OF  ENERGY 
DOE  PROGRAM  MANAGERS:  M.  MCCABE  AND  D.  ABRAMSON 
JANUARY  23.  1989 

COMMON/PRINTD/IDIA, IDIA0.N.M 

COMMON/PR I NT/ 1 PR 

COMMON/RESTR/IEXP 

COMMON/CONST/TC . PC . VC . TFR , AJ . EEP 

CONftADN/TGPG/AG . BG , CG . DG . EG . FG , AA . BB 

COMMON/DENSF/AL.BL.CL,DL,EL,BPL.CPL,DPL. EPL 

COMMON/ST AT E/A 1 ,B1 ,C1 . A2 . B2 , C2 , A3 , B3 .C3 , A4 . B4 , C4 . A5 , B5 , C5 . 

$ A6,B6,C6. ALPHA. AK 
COMMON/SPHTV/AC . BC . CC , DC . EC . FC . XX . YY 
C0f^0N/C0EFPR/A(5. 12) 

COMMON/HPHX/NSLABS.NDEP(2) .NROW(2) .DI (2) .DO(2) .DT(2) .TPCH(2) . 
k DPCH(2) .WIDTH(2) .FPCH(2) .FTK(2) .FMK(2) .TMK(2) .CFM1 .BSIDE(2) . 
k NTUB(2.5) . IFROM^.  130)  .NTPS(2)  .BSPACE(2) . 
k ACFM(2) . IFIN(2) . ISUR(2) .SFLOW{2) 

COMMON/AT EST/X ( 2 . 1 8 ) . VX U . 1 8 ) . NT EST ( 2 ) 

COt^ON/MERG/MERGE(2.20.2) . IMER(2) . IOUT(2 .20) .NOUT(2) . 
tIDEPTH(2. 130) .FLOW(2. 130) .KFEED(2 . 1 30 . 3) . KSTART (2 . 1 30) .KST(2) 
COMMON/MASS/PRM (2.2.130) 

COMMON/A I RD/I GET (2. 130.2) .GET(2. 130.2) .ARATIO(2. 130) 
COMMON/OUTPR/HR (2.130) 

DIMENSION  DATE(10) .RMS( 10) .HDS(10) .P2S( 10) .OLS( 10) 

OPEN(UNIT=5.  FI LE=’ INPUT') 

OPEN(UNIT=6.  F I LE=’ OUTPUT ■ ) 

OPEN(UNIT=10.  FILE=’DIAG’.  STATUS= 'UNKNOWN ’ ) 


WRITE(* .•) 'DATE: ' 
READ(5.800)WW1 .WW2.WW3.WW4 


800  FORMAT (4A4) 

WRITE(* .*) 
WRITE(*.*)’IPR  = 0 

WRITE(*.*)'IPR  = 1 
WRITE(* .•) ’ IDIA=  0 
WRITE(» .*) ' IDIA=  1 
1 ■ TOGETHER  WITH  MAIN 


A DEFAULT  DEVICE' 
TO  FILE  "RESULT"' 


FOR  MAIN  RESULTS  OUTPUT  ON 
FOR  MAIN  RESULTS  OUTPUT  IN 
FOR  NO  ADDITIONAL.  DETAILED  RESULTS  OUTPUT' 
FOR  ADDITIONAL.  DETAILED  RESULTS  OUTPUT'. 
RESULTS' 


WRITE(* . •) ’ IDIA=  2 FOR  ADDITIONAL.  DETAILED  RESULTS  OUTPUT 
1 ' FILES  "DIAG"  AND  D1AG0"' 

WRITE(» . •) ■ IPR.  IDIA  =■ 

READ(5.»)IPR. IDIA 
IF(IPR.EO.0)THEN 
1PR=6 
ELSE 


TO' 


IPR=9 

OPEN(UNIT=9.  FI LE=' RESULT ' . STATUS= ' UNKNOWN ' ) 
END  IF 

IF(IDIA.EQ.1)THEN 
IDIA=IPR 
ID1A0=IPR 
END  IF 

IF(ID1A.E0.2)THEN 

ID1A=10 


ID1A0=1 1 

OPEN  (UNIT=11.  FILE='DIAG0' . STATUS= ' UNKNOWN ' ) 

END  IF 
WRITE(* .*) 

WRITE(* . •) 'TAIR  = AIR  DRY  BULB  TEMPERATURE.  (F)' 

WRITE(* .•) 'RH  = AIR  RELATIVE  HUMIDITY.  (DECIMAL  FRACTION)' 
WRITE(* . •) 'TAIR.  RH  =' 

READ(5.*)TAIR.RH 
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WRlTEf*.*) 

WRITEf*.*) 'Xl  = REFRIGERANT  INLET  QUALITY,  (DECIMAL  FRACTION)’ 
WRITE(*.*)’TSAT2=  REFRIGERANT  SAT.  TEMP.  AT  COIL  OUTLET,  (F)’ 
WRITE(*,*) ■TSUP2=  REFRIGERANT  SUPERHEAT  AT  THE  COIL  OUTLET,  (F)’ 
WRITE(*,*) ’XI , TSAT2,  TSUP2  =’ 

READ(5,*)X1 ,TSAT2,TSUP2 


EVSIM  (VER.  1.1);  SIMULATION  OF  AN  EVAPORATOR’ 


WRITE(IPR,*) 

WRITE(IPR,900) 

WRITE(IPR,  •)’ 

1’  COIL’ 

WRITE(IPR,900) 

WRITE(IPR,801)WW1 ,WW2,WW3,WW4 
801  FORMAT (1X,4A4) 

C..*.  input  EVAPORATOR  COIL  DATA 

WRITE(IPR,*)’COIL  INFORMATION:’ 

CALL  RDATA3 

IF(IDIA.NE.0)  WRITE(IDIA0,*) ’•READING  OF  DTEV  COMPLETED’ 
WRITE(IPR,860)NSLABS 

860  FORMAT (’  NUMBERS  OF  SLABS  IN  THE  ASSEMBLY 1 14) 

WRITE(IPR,862)IEXP 

862  FORMAT(’  NUMBER  OF  EXPANSION  DEVICES: 

NCFM=CFM1 

WRITE(IPR,864)NCFM 
864  FORMAT (’  TOTAL  CFM: 

WRITEnPR,*)’AIR  CONDITION;’ 

WRITE(IPR,866)TAIR 
WRITE(IPR,868)RH 

866  FORMAT(’  AIR  DRY  BULB  TEMPERATURE: 

868  FORMAT^  AIR  RELATIVE  HUMIDITY: 

C.***  input  REFRIG.  PARAMETERS;  WRITE  THE  NAME  OF  REFRIGERANT  USED. 
CALL  DATAIN 

IF(IDIA.NE.0)  WRITE(IDIA0.*) ’•RETURNED  FROM  DATAIN’ 
WRITEnPR,^)  ’REFRIGERANT  CONDITIONS;  ’ 

WRITEnPR,870)X1 
WRITE! IPR.872)TSAT2 
WRITE(IPR,874nSUP2 

FORMAT!’  REFRIG.  QUALITY  AT  INLET; 

FORMATS  REFRIG 


MI4) 


’,116) 


’ . 1 F7 . 2 , 
• .ir7.2) 


F‘) 


870 

872 

874 


• .1F7.2) 

SAT.  TEMPERATURE  AT  EXIT : ’ . 1 F7 . 2 , ’ F’) 
FORMAT(’  REFRIG.  SUPERHEAT  AT  EXIT:  •.1F7.2.’  F’) 

DETERMINE  SIMULATION  CONTROLING  PARAMETERS 
IF(NSLABS.EQ.1)THEN 
SFLOW(1)=1 . 

IEXP=1 

KSLABS=1 

ELSE 

IF(IEXP.EO. 1)THEN 
SFLOW!1)=0.5 
SFLOW(2)-0.5 
KSLABS=2 
ELSE 

SFLOW(1)=1 ./(I .+SFLOW(2)/SFLOW(1)) 

SFLOW(2)=1 .-SFLOW(1 ) 

KSLABS=1 
END  IF 
END  IF 

DETERMINE  REFRIGERANT  CIRCUITRY  AND  DISTRIBUTION 
CALL  TRACE3 

IF(IDIA.NE.0)  WRITE(IDIA0,^) ’•RETURNED  FROM  TRACE3’ 

DETERMINE  AIR  DISTRIBUTION 
CALL  DISTR2 

IF(IDIA.NE.0)  WRITE(IDIA0,^) ’•RETURNED  FROM  DISTR2’ 


••••  ESTIMATE  REFRIGERANT  MASS  FLOW  RATE,  RMASS 
AREA=3.14^DI(1)^^2/4. 

G=1 .5E+5 

RMASS=G^AREA^REAL(NOUT( 1 ) ) •REAL(NSLABS) 
RMASS=AMAX 1 ! RMASS ,150.) 

RMASS=AM I N 1 ( RMASS , 680 . ) 

CALC.  ENTERING  REFRIG.  PRESSURE  AND  ENTHALPY 
T1=TSAT2+3. 

CALL  VPVHS(1  ,T1  ,P1 ,V1 .HI .SI .HF1 ) 
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H1=HF1+X1*(H1-HF1) 

C....  CALC.  REQUIRED  EXIT  ENTHALPY.  H2RE0 
PSAT=SATP(TSAT2) 

T2REO=TSAT2+TSUP2 

CALL  VPVHS(2 . T2REQ . PSAT , VG2 . H2RE0 . SG2 . HFG) 
CALL  CPCV( T2RE0 . PSAT . CP . CV , GA . SO) 

PSLOPE=1 . 


START  ITERATIONS 

IF(IDIA.NE.0)  WRITE(IDIA0.*)’*DATA  PREPARATION  COMPLETED’ 

C*  • • * 

IF(IDIA.NE.0)  WRITE(IDIA0.*) ’••••ITERATION  STARTS^^^^’ 

DO  110  N=1 , 10 
DO  107  110=1 , lEXP 
DO  105  M=1 ,5 

IF(IDIA.NE.0)  WRITE(IDIA0,^) ’N=  ’ .N. ’ M=’ .M. ’ SLAB#  ’,110 


CALL  EVPHX2(KSLABS. II0.RMASS.T1 .PI . HI . TAIR , 14 . 7 . RH . 
& T2.P2.H2,X2.S2.0L) 


TG2=SATT(P2) 

IF(IDIA.NE.0)  WRITE(IDIA.^) ’^702  =’.TG2 
PD=P2-PSAT 

IF(ABS(TG2-TSAT2).LT.0.05)GOTO  106 
IF(M.EQ. 1)THEN 

P1N=P1-PD^PSL0PE 

ELSE 

PSL0PE=(P1A-P1 )/(PDA-PD) 

pin=pi-pd^pslope 

END  IF 
P1A=P1 
P1=P1N 
T1=SATT(P1) 

105  PDA=PD 

WRITE(IPR.^) ’ loop  105  DID  not  CONVERGE,  HO. PD  =’,110, PD 

106  IF(IEXP.EQ.2)THEN 

IF(I10.EQ.1)THEN 

H22=H2 

P22=P2 

0L1-QL 

ELSE 

H2=H22^SFL0W(1)+H2^SFL0W(2) 

P2=P22^SFL0W(1 )+P2^SFLOW(2) 

0L=0L+QL1 
END  IF 
END  IF 

107  CONTINUE 
H2DIF=H2-H2REQ 

1F(IDIA.NE.0)  WRITE(IDIA.^)’^^^^^^^^^H2DIF=’ ,H2DIF 

RMS(N)=RMASS 

HDS(N)=H2DIF 

OLS(N)=QL 

P2S(N)=P2 

IF(N.GT.2)THEN 

I F ( ABS ( HDS ( N- 1 ) ) . LT . 5 . • CP . AND . ABS ( H2D I F ) . LT . 2 . •CP ) GOTO  1 50 
END  IF 

1F(N.EQ.1.)THEN 

RMASSN=RMASS^ (H2-H1 )/(H2RE0-H1 ) 

ELSE 

RMASSN=RMASSA-H2D I FA • ( RMASS-RMASSA ) / ( H2D I F-H2D I FA ) 

END  IF 

IFfRMASSN.GT.RMASS)RMASSN=AMIN1 (RMASSN. 1 .5^RMASS) 

I F ( RMASSN . LT . RMASS ) RMASSN=AMAX 1 ( RMASSN . 0 . 7 • RMASS ) 

P 1 -PSAT+ ( P 1 -P2 ) • ( RMASSN/RMASS ) • • 2 
RMASSA=RMASS 
RMASS=RMASSN 
110  H2DIFA=H2DIF 

IF(IDIA.NE.0)  WRITEdDIA,  •)'••• 

WRITE(IPR.^) ’LOOP  110  DID  NOT  CONVERGE’ 

C^^^^  CALC,  k PRINT  RESULTS  FOR  THE  LAST  ITERATION  LOOP 
150  0T=RMASS^(H2-H1) 

DTIT=2. 

PP2-P2 
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CALL  PHWET2(PP2.H2.DTIT. 0.001 , T2 , V2 . S2 . X2 . TG2) 

TSUP=T2-TG2 

WRITE(IPR.900) 

WRITEf IPR,*) 

WRITE(IPR.  •)'  SIMULATION  RESULTS’ 

WRITE(IPR.*) 

WRITE(IPR.900) 

900  FORMAT{ 

WRITEnPR,*)’ •RESULTS  FROM  THE  LAST  ITERATION  LOOP:’ 

WRITEflPR.*) 

WRITE(IPR.905)QT 

905  FORMAT(’  TOTAL  CAPACITY:  ’,1F7.0.’  BTU/H’) 

WRITE(IPR.907)OL 

907  FORMAT (’  LATENT  CAPACITY: ’. 1 F7 .0, ’ BTU/H’) 

WRITE(IPR.910)RMASS 

910  FORMAT(’  REFRIGERANT  MASS  FLOW  RATE : ’ . 1 F6 . 1 . ’ LB/H’) 

WRITE(IPR.*) ’ REFRIGERANT  PARAMETERS  AT  THE  COIL  OUTLET:’ 
WRITE(IPR.915)TG2 

915  FORMAT (9X. ’SATURATION  TEMPERATURE :’. 1 F6 . 1 . ’ F’) 
WRITE(IPR.920)TSUP 

920  FORMAT(22X, ’SUPERHEAT: ’ .1F6.1 . ’ F’) 

WRITE(IPR.925)P2 

925  FORMAT(23X. ’PRESSURE: ’ .1F6. 1 . ’ PSIA’) 

WRITE(IPR,930)H2 

930  FORMAT (23X. ’ENTHALPY: ’ ,1F6.1 , ’ BTU/LB’) 

C****  PREPARE  AND  PRINT  REFRIGERANT  INFO  FOR  INDIVIDUAL  OUTLET  TUBES 
WRITE(IPR.*) 

WRITE! IPR,*) ’ INFORMATION  ON  REFRIGERANT  LEAVING  OUTLET  TUBES’ 
WRITE! IPR,*) 

WRITE(IPR,940) 

940  FORMAT ( ’ SLAB  # ’ . 3X . ’ TUBE  # ’ . 6X . ’ P ’ , 8X . ’ T ’ . 6X . ’ TSUP ’ , 6X . ’ X ’ . 7X , 

1 ’RMS’,/.’  (-) ’ ,6X. ’ (-) ’ .5X. ’ (PSIA) ’ .5X. ’ (F) ’ .6X. ’ (F) ’ .5X. 

2 ’(-)’.4X.’(LB/H)’) 

DO  160  IIO=1.NSLABS 
DO  160  L=1 .NOUT(IIO) 

I=IOUT(IIO,L) 

HOUT=HR(IIO. I) 

POUT=PRM(IIO,2,I) 

CALL  PHWET2 ( POUT. HOUT.DT IT. 0.001 .T2.V2.S2.X2.TG2) 
RMSI=RMASS*FLOW(IIO. I) 

TSUP=T2-TG2 

160  WRITE(IPR.944)II0. I .POUT.T2.TSUP.X2.RMSI 

944  FORMAT(I5. IX. 19. 1X.1F10.2.1F8.1 . 1 F9 . 1 . 1 F9 . 3 . 1 F9 . 2) 

C.*.*  CALC.  AND  PRINT  RESULTS  FOR  REQUESTED  SUPERHEAT 
IF(H2DIF.NE.0. )THEN 
HDH=HDS(N) 

RMSH=RMS(N) 

P2H=P2S(N) 

QLH=QLS(N) 

HDL=HDS(N-1 ) 

RMSL=RMS(N-1 ) 

P2L=P2S(N-1) 

QLL-QLS(N-I) 

DO  120  1=1 .N-2 
IF(HDH.NE.HDL)GOTO  122 
HDL=HDS(N-1-I) 

RMSL=RMS(N-1-I) 

P2L=P2S(N-1-n 

120  0LL=QLS(N-1-I) 

122  H2L=H2RECH-HDL 

H2H=H2RECH-HDH 

HS  LOPE= ( H2R  EO-H2  L ) / ( H2H-H2  L ) 

RMASS=RMS  L+ ( RMSH-RMS  L ) • HS  LOP  E 
P2=P2  L+ ( P2H-P2  L ) • HS  LOP  E 
QL=OLL+(QLH-QLL)*HSLOPE 
H2=H2REQ 
END  IF 

QT=RMASS*(H2-H1 ) 

DTIT=2. 

PP2=P2 

CALL  PHWET2(PP2.H2.DTIT. 0.001 .T2.V2.S2.X2.TG2) 

TSUP=T2-TG2 


74 


c*...  PRINT  RESULTS 
WRITE(IPR.*) 

writeopR.900) 

WRITE(IPR,9e2hsUP2 

902  FORMATC  •RESULTS  FOR  THE  REQUESTED  REFRIGERANT  SUPERHEAT F5 . 1 ) 
WRITE(IPR,*) ’ (INTERPOLATED  FROM  RESULTS  OF  LAST  TWO  ITERATIONS)’ 
WRITE(IPR,») 

WRITE(IPR,905)OT 
WRITE(IPR,907)OL 
WRITE (I PR  910)RMASS 

WRITE(IPR,’*)  ’ REFRIGERANT  PARAMETERS  AT  THE  COIL  OUTLET:’ 

WRITE(IPR,915)TG2 

WRITE(IPR.920)tSUP 

WRITE(IPR.925)P2 

WRITEf IPR,930)H2 

WRITE(IPR,900) 

STOP 

END 
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SUBROUTINE  AIRPR(I .T.PATM.RH.W.CP.R.AM.AK) 

•••*  PURPOSE: 

TO  CALCULATE  AIR  PROPERTIES 
JANUARY  17,  1989 

INPUT  DATA: 

I = 1 IF  RELATIVE  HUMIDITY  IS  GIVEN  (-) 

=2  IF  HUMIDITY  RATIO  IS  GIVEN  (-) 

T - AIR  TEMPERATURE  (-) 

PATH  - AIR  PRESSURE  (PSIA) 

RH  - RELATIVE  HUMIDITY  IN  FRACTION  (IF  1=1)  (-) 

W - HUMIDITY  RATIO  (IF  1=2)  (LBM  H20/LBM  DRY  AIR) 

•••*  OUTPUT  DATA: 

AK  - AIR  THERMAL  CONDUCTIVITY  (BTU/H*F*FT) 

AM  - AIR  DYNAMIC  VISCOSITY  (LBM/FT*H) 

CP  - AIR  SPEC.  HEAT  AT  CONSTANT  PRESSURE  (BTU/LBM*F) 
R - GAS  CONSTANT  OF  AIR  ( LBF*FT/LBM*R) 

RH  - RELATIVE  HUMIDITY  IN  FRACTION  (FOR  1=2)  (-) 

W - HUMIDITY  RATIO  (FOR  1=1)  (LBM  H20/LBM  DRY  AIR) 


DOUBLE  PRECISION  TR.Z.PSAT.P.WD.RHD.PW.CPD 
TR=T+460 . 

Z=1000./TR 

IF(TR.GE.492.)GOTO10 

PSAT=DEXP(0 . 03940*Z**3-0 . 2755*Z*Z-1 0 . 431  •Z+l  9.509) 

GOTO30 

10  IF(TR.GE.672.)GOTO20 

PSAT=DEXP(0 . 1 7829*Z*Z*Z-1 . 6896*Z*Z-5 . 0988*Z+1 3 . 4353) 

GOTO30 

20  PSAT=DEXP(0.71692*Z**4-4.01506*Z**3+7.5568*Z*Z-14.2131*Z+16.8255) 
30  IF(I .EQ.2)THEN 

P=W* PATM/ ( 0 . 6 1 98D0+W ) 

RHD=P/PSAT 

RH-RHD 

WD-W 

ELSE 

RHD*RH 

I F ( RH . GE . 0 . 0000 1 ) GOTO40 
W=0. 

WD=0 . 00 
GOTO50 

40  PW=RHD*PSAT 

WD-0 . 61 98*PW/(PATM-PW) 

W=WD 

50  CONTINUE 

END  IF 

CPD=0 . 2478786-0 . 4204563E-04*TR+0 . 5767857E-07*TR**2 
k -0. 1493056E-10*TR**3 
CP-(CPD+0.444*WD)/(1 .+WD) 

R=(53. 34+85. 76*WD)/(1 .4WD) 

AM=5 . 5029E-03+8 . 71 57E-05*TR-2 . 9464E-08*TR**2 
k +6.250E-12*TR**3 

AK— 2 . 853E-04+3 . 268E-05*TR-8 . 253E-09*TR*TR 
k +1 .239E-12*TR**3 
RETURN 
END 
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FUNCT I ON  A I RHT3 ( I ROW , DO , TPCH . DPCH . FPCH , FTK , W I DTH , DH . NDEP . I F I N . 
k AMASS. AVIS. ACP.AK) 

••••  PURPOSE:  TO  CALCULATE  DRY  AIR-SIDE  HEAT  TRANSFER  COEFFICIENT 
FOR  A TUBE  WITH  FLAT.  WAVY.  AND  LANCED  PLATE  FINS. 
SEPTEMBER  22  .1988. 

••••  NOTE:  THIS  FUNCTION  CALCULATES  H.T.C  FOR  A TUBE  DEPENDING 
ON  THE  DEPTH  LOCATION  OF  THE  TUBE  IN  THE  SLAB. 

GEOMETRIES  ARE  ASSUMED  FOR  WAVY  k STRIP  FINS. 

••••  INPUT  DATA: 

ACP  - AIR  SPEC.  HEAT  AT  CONSTANT  PRESSURE  (BTU/LBM*F) 

AK  - AIR  THERMAL  CONDUCTIVITY  (BTU/H*FT*F) 

AMASS  - AIR  MASS  FLOW  RATE  ASSOCIATED  WITH  THE  TUBE  (LBM/H) 

AVIS  - AIR  DYNAMIC  VISCOSITY  (LBM/FT*H) 

DH  - HYDRAULIC  DIAMETER  OF  THE  SLAB  (FT) 

DO  - TUBE  OUTSIDE  DIAMETER  (FT) 

DPCH  - TUBE  DEPTH  PITCH  (FT) 

FPCH  - FIN  PITCH  (FT) 

FTK  - FIN  THICKNESS  (FT) 

I FIN  - FIN  DESIGN  DESCRIPTOR  (-) 

■=1  FOR  A FLAT  FIN 
-=2  FOR  A WAVY  FIN 

-=  3 FOR  AN  ENHANCED  FIN  (LANCED.  STRIPPED) 

IROW  - DEPTH  ROW  LOCATION  OF  THE  TUBE  (-) 

NDEP  - NUMBER  OF  DEPTH  ROWS  IN  THE  SLAB  (-) 

TPCH  - TUBE  PITCH  IN  THE  ROW  (FT) 

WIDTH  - HEAT  EXCHANGER  WIDTH  (TUBE  LENGHT)  (FT) 

•••  OUTPUT  DATA: 

AIRHT3  - DRY  AIR-SIDE  HEAT  TRANSFER  COEFFICIENT  (BTU/H*F*FT**2) 

COMMON/PR I NT/ I PR 
REAL  JM. JN,J4. JTUBE.NUF.NUW 
DATA  PI/3. 1415927/. IERR/0/ 

•••  FLAT  FIN  •••*••***•»*••*•••**•**••••••**•**••••**•*•*••••••••• 

•••  REFERENCE:  GRAY.  D.L.  AND  WEBB.  R.L..  HEAT  TRANSFER  AND  FRICTION 
CORRELATIONS  FOR  PLATE  FINNED-TUBE  HEAT  EXCHANGERS 
HAVING  PLAIN  FINS.  PROC.  OF  EIGHTH  INT.  H.T. 

CONFERENCE.  SAN  FRANCISCO.  1986. 

FN-WIDTH/FPCH 

AREAC-(TPCH-DO) • (WIDTH-FN*FTK) 

GOAMASS/AREAC 
RED-GC*DO/AVIS 

I F (RED . LT . 500 . . OR . RED . GT . 24700 . )THEN 
WRITE(10.*)*AIRHT3.  RED  LIMIT  500-24700.  RED  ='.RED 
END  IF 

SPACE-FPCH-FTK 
TR-REAL(IROW) 

TR-=AMIN1(TR.4.) 

J4-0.14*RED**(-. 328)* (TPCH/DPCH)**(-. 502)* (SPACE/DO)**0. 0312 
JN-:J4*0.991*(2.24*RED**(-.092)*(TR/4.)**(-.031))**( .607*(4.-TR)) 
IF(TR.EQ.1 .)THEN 
JTUBE=JN 
ELSE 

JM=(2.24*RED**(-.092)*((TR-1 .)/4. )••(-. 031 ))••(. 607* (5. -TR)) 
jm=J4*0.991*JM 
JTUBE<=TR*JN-(TR-1 . )*JM 
END  IF 

PR»AVIS*ACP/AK 

AIRHT3«JTUBE*GC*ACP/PR**. 66667 
IF(IFIN.E0.3)THEN 

• ••  LOUVERED/LANCED  FINS  

•••  REFERENCE:  NAKAYAMA.  W.  AND  XU.  L.P. . ENHANCED  FINS  FOR  AIR- 
COOLED HEAT  EXCHANGERS  - HEAT  TRANSFER  AND  FRICTION 
FACTOR  CORRELATIONS.  ASME-JSME  THERMAL  ENGINEERING 
JOINT  CONFERENCE.  PROCEEDINGS.  PP.  495-510.  ASME,  NY. 
MARCH  1983. 

REH-GC*DH/AVIS 

C****  NOTE:  CORRELATION  IS  GOOD  FOR  FIS  IN  THE  RANGE  0.2  - 0.35 
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FIS=0.275 

IF(IERR.EQ.0)THEN 

IERR=1 

GAP=SPACE*3048. 

IF(GAP. LT.0.15.OR.GAP.GT.0.2)WRITE(10,*) •AIRHT3. 

1 GAP  LIMIT  IS  0.15  - 0.2  MM  . GAP  «=  '.GAP 
FP=FPCH*3048. 

IF(FP. LT. 1 .8.OR.FP.GT.2.5)WRITE(10,*) ■AIRHT3. 

1 FIN  PITCH  LIMIT  IS  1.8  - 2.5  m.  FP  «='.FP 

IF(REH.LT.250. .OR.REH.GT.3000. )WRITE(10.*) 'AIRHTS. 

1 REH  LIMIT  IS  250  - 3000,  REH  =',REH 
END  IF 

FJ-=1 .+1 093 . * ( FTK/SPACE) •• 1 . 24*FIS**0 . 944*REH** (-0 . 58)+ 

1 1 .097*(FTK/SPACE)**2.09*FIS**2.26*REH**0.88 

AIRHT3-=AIRHT3*FJ 
END  IF 

IF(IFIN.EQ.2)THEN 

WAVY  FINS  

••••  REFERENCE:  TRAUGER,  P.  AND  WEBB.  R.L..  A CORRELATION  FOR  THE 
AIR  SIDE  HEAT  TRANSFER  COEFFICIENT  FOR  A WAVY  FIN. 
TO  BE  PUBLISHED,  1988. 

THIS  CORRELATION  WAS  DEVELOPED  FOR  3 DEPTH  ROW  H-X. 
••••  ASSUMPTIONS:  1.  TWO  FIN  WAVES  PER  DEPTH  PITCH 

2.  WAVE  DEPTH  EQUALS  THE  SPACE  BETWEEN  FINS 
DC>=DC>+2.*FTK 
SP2=DPCH/2 . 

SD=SPACE 

BETA-=0 . 25*PI  •DC*DC/(TPCH*DPCH) 

GAMMA=SQRT(1 .+4. •(SD/SP2)**2) 

DHW=2. -SPACE*  (1  .-BETA)/(GAIwMA*(1  .-BETA)+2 . -SPACE-BET  A/DC) 
DHF-=2. -SPACE* (1  .-BETA)/(  1 .-BETA+2 . -SPACE-BETA/DC) 
REDF-=AMASS*DHF/(AVIS*AREAC*  ( 1 . -BETA)  ) 

R EDW=AMASS * DHW/ ( A V I S * AR EAC * ( 1 . -B ET A ) ) 

DEP=REAL(NDEP) 

GZW=REDW*PR*DHW/fDPCH*DEP) 

GZF-=REDF*PR*DHF/(DPCH*DEP) 

IF(GZW.LT.25.)THEN 

NUW«0 . 5-GZW*  * . 86* ( TPCH/DC ) * * . 1 1 * ( SPACE/DC )**(-. 09 ) 
NUW=NUW* (SD/DPCH) ** . 1 2* (SP2/DPCH) * * (- . 34) 

NUF-0 . 4-GZF--0 . 73* (SPACE/DC) ** (-0 . 23) *3 . * *0 . 23 
ELSE 

NUW»0 . 83-GZW*  *0 . 76* ( TPCH/DC ) * * . 1 3* ( SPACE/DC )**(-.  1 6 ) 
NUW=NUW* (SD/DPCH) ** . 25* (SP2/DPCH) * * (- . 43) 

NUF-0 . 53-GZF*  *0 . 62* (SPACE/DC) ** (-0 . 23) *3 . **0 . 31 
END  IF 
CF-NUF/GZF 
CW-NUW/GZW 

IF(CF.GE. .5.0R.CW.GE. .5)THEN 

WRITE(10.*)'AIRHT3.  IFIN=2.  TOO  LOW  AIR  FLOW  RATE* 
WRITEO0.*) ’fLAT  FIN  H.T.C  ASSUMED* 

ELSE 

NUF-.25*GZF*ALOG((1 .+2 . *NUF/GZF)/( 1 .-2.-NUF/GZF)) 
NUW=.25*GZW*AL0G(O .+2 . *NUW/GZW)/( 1 .-2 . -NUW/GZW) ) 
FJ-NUW/NUF 
AIRHT3»AIRHT3*FJ 
END  IF 
END  IF 
RETURN 
END 
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SUBROUTINE  BALFL1 (KSLABS, 1 1 1 ) 


••••  PURPOSE: 

TO  ADJUST  REFRIGERANT  FLOW  DISTRIBUTION  IN  AN  EVAPORATOR 
BASED  ON  PRESSURE  DROP  AT  DIFFERENT  CIRCUITS. 

THE  PROGRAM  CAN  HANDLE  SINGLE  SLAB  AND  TWO  SLAB  (A  k V)  COILS. 
10-14-1988 


••••  INPUT  DATA: 

IFROM(IIO.J) 


III 

KSLABS 

IMER(IIO) 

lOUT(IIO.L) 

KFEED(IIO.J.N) 


KSTART(IIO.N) 

KST(IIO) 

MERGE(IIO.K.I) 

MERGE(IIO.K.2) 

NDEP(IIO) 

NOUT(IIO) 

NTPS(IIO) 
PRM(II0.1 .1) 
PRMUI0.2.I) 


NUMBER  OF  THE  TUBE  FROM  WHICH  TUBE  J RECEIVES 
REFRIGERANT.  IF  THE  TUBE  IS  CONNECTED  TO  THE 
INLET  MAINFOLD,  I FROM  IS  SET  TO  0. 

1 FOR  THE  FIRST  SLAB  (-) 

2 FOR  THE  SECOND  SLAB  (-) 

NUMBER  OF  SLABS  IN  THE  EVAPORATOR  ASSEMBLY 
TO  BE  CONSIDERED  (KSLABS-2  OVERRIDES 
SPECIFICATION  OF  III)  (-) 

NUMBER  OF  SPLIT  POINTS  (-) 

NUMBER  OF  THE  TUBE  CONNECTED  TO  THE  OUTLET 
MANIFOLD.  FOUND  AS  L SUCH  TUBE  (-) 

NUMBER  OF  THE  TUBE  RECEIVING  REFRIGERANT 
FROM  TUBE  J.  FOUND  AS  N SUCH  TUBE. 

NOTE  THAT  TUBE  J CAN  FEED  UP  TO  3 TUBES 
(N  CAN  BE  1.2  AND  3).  KFEED  IS  SET  TO  -1  IF 
J TUBE  FEEDS  THE  DISCHARGE  MANIFOLD.  KFEED 
IS  SET  TO  0 IF  A TUBE  IS  NOT  FED.  (-) 

NUMBER  OF  THE  TUBE  CONNECTED  TO  THE  INLET 
MANIFOLD,  FOUND  AS  N SUCH  TUBE  (-) 

NUMBER  OF  TUBES  CONNECTED  TO  THE  INLET 
MANIFOLD  (-) 

NUMBER  OF  THE  TUBE  WHICH  FEEDS  A SPLIT  POINT. 
FOUND  AS  K SUCH  TUBE  (-) 

NUMBER  OF  TUBES  FED  BY  THE  TUBE  K (-) 

NUMBER  OF  TUBE  DEPTH  ROWS  IN  THE  SLAB  (-) 
NUMBER  OF  TUBES  CONNECTED  TO  THE  OUTLET 
MANIFOLD  (-) 

NUMBER  OF  TUBES  IN  THE  SLAB  (-) 

REFRIG.  PRESSURE  AT  INLET  OF  TUBE  I (PSIA) 
REFRIG.  PRESSURE  AT  OUTLET  OF  TUBE  I (PSIA) 


•••  OUTPUT  DATA: 

FLOW(IIO.J)  - FRACTION  OF  COIL  TOTAL  REFRIG.  MASS  FLOW 
PASSING  THROUGH  TUBE  J (-) 


•••  SUBPROGRAMS  CALLED  BY  BALFL1 : FRACT 

COMMON/HPHX/NSLABS.NDEP(2) .NROW(2) ,DI(2) ,DO(2) .DT(2) .RPCH(2) . 
k DPCH(2) .WIDTH(2) .FPCH(2) .FTK(2) .FWB<(2) .TMK(2) .CFM1 ,BSIDE(2) . 
k NTUB(2.5) , IFROM(2, 130) ,NTPS(2) .BSPACE(2) . 
k ACFM(2) . IFIN(2) . ISUR(2) ,SFLOW(2) 

COMwlON/MERG/MERGE(2.20.2)  . IMER(2)  . IOUT(2.20)  .NOUT(2)  . 
k IDEPTH(2. 130)  .FLOW(2. 130KKFEEDU.  130,3), KSTART(2. 130). KST(2) 
COH(WON/MASS/PRM  (2,2.130) 

DIMENSION  LEFT(20) .PC(2, 130) ,RN(20) , F(20) . ITUBE(20) , ISEE(20) 


•••  FIND  REFRIGERANT  FLOW  DISTRIBUTION 


Ill-Ill 

DO  120  IDO=1, KSLABS 
IF(KSLABS.EQ.1)THEN 
IIO-II1 
ELSE 

IIO-IDO 
END  IF 

DO  65  1-1 .IMER(IIO) 

65  LEFT(I)-MERGE(IIO.I.2) 
DO  70  1-1 .NTPS(IIO) 

70  PC(IIO.I)-0. 

DO  120  IL-1 .NOUT(IIO) 

I-IOUT(IIO.IL) 

P0UT-PRM(II0.2.I) 

DO  100  IT-1 .NTPS(IIO) 
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IP-=IFROM(nO.  I) 

IF(IP.EQ.C)THEN 
PC(IIO. I)-POUT 
GOTO  120 
END  IF 

IF(KFEED(IIO.IP.2).EQ.0)THEN 

I=IP 

GOTO  100 
END  IF 

PC(IIO.I)-POUT 
DO  75  IM=1 .IMER(IIO) 

75  IF(IP.EQ.MERGE(IIO.IM.1))GOTO  77 
77  LEFT(IM)-LEFT(IM)-1 

IF(LEFT(IM).GT.0)GOTO  120 
POUT-10. 

PUP=PRM(II0.2,IP) 

NSPLIT-MERGE(IIO.IM.2) 

DO  90  11=1 .NSPLIT 
N=KFEED(IIO.IP.I1) 

POUT«POUT+PC(IIO.N)*FLOW(IIO.N)/FLOVy(IIO.IP) 

90  RN(I1)-=(PUP-PC(IIO.N))/FLOW(IIO.N)**1 .75 
CALL  FRACT(NSPLIT.RN,F) 

DO  92  11=1 .NSPLIT 
N=KFEED(IIO.IP.I1) 

92  FLOW(IIO.N)-F(I1) 

I-IP 

100  CONTINUE 
120  CONTINUE 

IF(KSLABS.EQ.1)THEN 

NSTART=KST(IIO) 

ELSE 

NSTART-KST ( 1 )+KST ( 2 ) 

END  IF 

DO  130  I=1.NSTART 
IF(KSLABS.EQ.1)THEN 
IIO-II1 

N=KSTART(IIO.I) 

ELSE 

IF(I.LE.KST(1))THEN 

110=1 

N=KSTART(1 .1) 

ELSE 

110=2 

N1-KST(1) 

N-KSTART(2.I-N1) 

END  IF 
END  IF 

130  RN(I)=(PRM(II0.1 .N)-PC(II0,N))/FL0W(II0.N)**1 .75 
CALL  FRACT(NSTART,RN,F) 

DO  132  I=1,NSTART 
IF(KSLABS.EQ.1)THEN 
N-KSTART(IIO.I) 

ELSE 

IF(I.LE.KST(1))THEN 

110=1 

N=KSTART(1 ,1) 

ELSE 

110=2 

N1-KST(1) 

N=KSTART(2,I-N1) 

END  IF 
END  IF 

132  FLOW(IIO.N)=F(I) 

*•••  ASSIGN  REFRIGERANT  DISTRIBUTION.  FLOW(IIO.I) 
ISTORE=0 
110=111 

DO  170  ID0=1.KSLABS 
IF(KSLABS.EQ.2)IIO=IDO 
DO  136  1-1 ,IMER(IIO) 

ITUBE(I)-0 
136  ISEE(I)-0 

DO  160  IS-1 .KST(IIO) 
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I-:KSTART(IIO.  IS) 

IL-1 

DO  150  10-1 .NOUT(IIO) 

DO  145  IT-1 ,NTPS(IIO) 

INI-KFEED(IIO.I.IL) 

IF(IN1 .EQ.-1)THEN 
IF(ISTORE.GT.0)THEN 
I-ITUBE(ISTORE) 

IL-ISEE(ISTORE) 

I STORE- 1 STORE- 1 
GOTO  150 
END  IF 
GOTO  160 
END  IF 

IF(IL.GT.1)G0T0  137 
DO  135  11-2,3 
IN2-KFEED(II0.I.I1) 

IF(IN2.EQ.0)GOTO  137 
ISTORE-ISTORE+1 
ITUBE(ISTORE)-I 
135  ISEE(IST0RE)-I1 
137  IN2-KFEED(IIO.I.2) 

IF(IN2.GT.0)THEN 

FLOW(IIO.IN1)-FLOW(IIO.IN1)*FLOW(IIO.I) 

ELSE 

FLOW(IIO. INf)»FLOW(IIO.I) 

END  IF 

I-IN1 

IL-1 

145  CONTINUE 
150  CONTINUE 
160  CONTINUE 
170  CONTINUE 

IF(KSLABS.EQ.1)THEN 
DO  180  1-1 .NTPS(IIO) 

180  FLOW(IIO.I)-FLOW(IIO.I)*SFLOW(IIO) 

END  IF 
RETURN 
END 
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SUBROUT I NE  CPCV ( TS . P . CP . CV . GAMMA . SON I C ) 

•••*  PURPOSE: 

TO  CALCULATE  FOR  REFRIGERANT  VAPOR 
SPECIFIC  HEAT  AT  CONSTANT  PRESSURE. 

SPECIFIC  HEAT  AT  CONSTANT  VOLUME. 

SPECIFIC  HEAT  RATIO  AND  SONIC  VELOCITY. 

****  JANUARY  9.  1989 

••••  INPUT  DATA: 

TS  - REFRIG.  VAPOR  TEMPERATURE  (F) 

P - REFRIG.  VAPOR  PRESSURE  (PSIA) 

REFRIG.  CONSTANTS  - REFER  TO  THE  MAIN  PROGRAM. 

••••  OUTPUT  DATA: 

CP  - SPEC.  HEAT  AT  CONSTANT  PRESSURE  (BTU/LBM*F) 

CV  - SPEC.  HEAT  AT  CONSTANT  VOLUME  (BTU/LBM*F) 

GAVMA  - SPEC.  HEAT  RATIO  (-) 

SONIC  - SONIC  VELOCITY  (FT/SEC) 

••••  SUBPROGRAMS  CALLED  BY  CPCV: 

SATT.VPSV 

DOUBLE  PRECISION  Z.Z2.Z3.V2.V3.V4,V5.V6.AKTTC.AXTTC.FDPDV. 

1 FDPDT.DPDT.FCCV.CVD.CPD.GAMMAD 
COMMON/PR I NT/ I PR 
COlwMON/CONST/TC . PC . VC . TFR . AJ  . EEP 

COMWION/STATE/AI  .B1  .Cl  .A2.B2.C2.A3.B3.C3,A4.B4.C4.A5.B5.C5. 
4:A6.B6,C6.ALPHA.AK 
COMMON/SPHTV/AC . BC . CC . DC . EC . FC . X . Y 
SAVE  TS  LAST . P LAST . CP  LAST , CVLAST . G LAST . S LAST 
DATA  TSLAST.PLAST/-1 . .-1 ./ 

TmTQj.TFR 

IFrT.LE.0.)GOTO  999 
IFfP.LT.e. )GOTO  999 
IF(ABS(TS-TSLAST).GT.  1.0E-4)GOTO  5 
IF(ABS(P-PLAST).GT.  1.0E-4)GOTO5 
CP-CPLAST 
CV-CVLAST 
GAMvIA-GLAST 
SONIOSLAST 
RETURN 
5 TG-SATT(P) 

IF(TS.LT.TG)TS«TG 

V- VPSV(P.TS) 

VI- V-B1 
V2-=V1  •VI 
V3«V1 •V2 
V4-V1 •V3 
V5=V1*V4 
V6-VUV5 
AKTTC»AK*T/TC 
AXTTC=DEXP(-AKTTC) 

Z=ALPHA*V 
Z2=2.*Z 
Z3=3.*Z 

IF(Z.GT. 150.D0)Z-150.D0 
IFrZ2.GT.150.D0)2-150.D0 
IF(Z3.GT. 150.D0)Z3-150.D0 
FDPDV=0 . 

FDPDT=0 . 

I F(ABS(A6) . LE. 1 . E-20)GOTO4 
IF(ABS(C1) .GE.1 .E-20)gOTO3 
FDPDV=-A  LPHA*  DEXP (-Z ) • ( A6+B6*T ) 

FDPDT=B6*DEXP(-Z) 

G0T04 

3 FDPDV=-(ALPHA*(DEXP(-Z3)+2.*C1*DEXP(-Z2))/(DEXP(-Z2)+2.*C1* 
acDEXP  (-Z  )+C  1 ‘C 1 ) ) • ( A6+B6*T+C6*  AXTTC  ) 

FDPDT=(B6-AK*C6*AXTTC/TC)*DEXP(-Z2)/(DEXP(-Z)+C1) 

4 DPDV-^AI  •T/V2-2 . • (A2+B2*T+C2*AXTTC)  A3-3 . • (A3+ 
«cB3*T4C3*AXTTC)  A++4  • • (A4+B4*T+C4*AXTTC)  A5-5  • • (A5+B5 
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a:*T4C5*AXTTC)A6+FDPDV 

DPDT-A 1 /V 1 + ( B2-AK  *02  • AXTTC/TC  ) A2+  ( B3-AK  *03  • AXTTC/TC  ) 
^^^(B4-AK•C4•AXTTC/TC)A^(B5-AK•C5*AXTTC/TC)A5+f^DPDT 
FCCV-0 . 

IF(ABS(C1 ) .GE. 1 .E-20)FCCV-C6*DEXP(-Z)/ALPHA-(C6*C1/ALPHA)* 
a:DL0G(1  .D0+DEXP(-Z)/C1) 

CVD-AC+BC • T+CC • T •• 2+DC • T •• 3+FC/T •* 2- ( 0 . 1 85053 • AK •• 2 • T •AXTTC/TC •• 2 ) 
k*  (C2Al+C3/(2.  •V2)+C4/(3 . •V3)4C5/(4.  •V4)+FCCV) 

CPD-CVD-0 . 1 85053*T*DPDT* •2/DPDV 

GAH(^tAD-CPD/CVD 

CP-CPD 

CV-CVD 

GAMMA-CAKMAD 

SONIOV*DSQRT(857 . 36091  •T*DPDT**2/CVD-4633 . 056*DPDV) 

CPLAST-CP 

CVLAST-CV 

GLAST-GAKA4A 

SLAST=S0NIC 

TSLAST-TS 

PLAST-P 

RETURN 

999  WRITE(IPR.100) 

100  FORMAT (5X. ’ERROR  IN  CALLING  -CPCV-’) 

RETURN 

END 
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SUBROUTINE  DATA IN 


•••*  PURPOSE: 

TO  READ  REFRIGERANT  CONSTANTS 
1/5/87 

COlyWON/PRINT/IPR 
COM^N/CONST/TC . PC . VC . TFR . AJ  . EEP 
COM^N/TGPG/AG , BG . CG . DG . EG . FG . AA , BB 
COt^N/DENSF/AL.BL.CL.DL.EL.BPL.CPL.DPL.EPL 
COMMON/ST AT E/A 1 .B1 .Cl .A2,B2,C2.A3.B3.C3,A4.B4.C4,A5.B5.C5. 
$ A6,B6.C6. ALPHA. AK 
COt«MON/SPHTV/AC . BC . CC . DC . EC . FC . X . Y 
COMMON/COE  FPR/A (5.12) 

OPEN  (UNIT=7.FILE-'DATAREF* . STATUS* ' 0 LD ' ) 



c****  INPUT  REFRIGERANT  DATA 

READ(7.800)T1 .T2.T3.T4.T5 
WRITE(IPR.802)T1 .T2.T3.T4.T5 
READ(7.*)TC.PC.VC 

read(7.*)tfr.aj.eep 

READf7.*)AG.BG.CG 

READ(7.*)DG.EG.FG 

READ(7.*)aA.BB 

READf7.*)AL.BL.CL 

READ(7.*)DL.EL.BPL 

READ{7..)CPL.DPL.EPL 

READ(7.*)a1 .B1 .Cl 

READ(7.*)a2.B2.C2 

READ(7.*)A3.B3.C3 

READ(7.*)A4.B4.C4 

READ(7.*)a5.B5.C5 

READ(7.*)A6.B6.C6 

read(7.*)alpha.ak 

read(7.*)ac.bc.cc 

READ(7.*)DC.EC.FC 

READ(7.*)X.Y 

DO10I-1 .5 

READ(7.*)fA(I .J) .J-1 .3) 

READ(7.*)(An.j).J-4.6) 

READ(7.*)(A( I .j) .J-7.9) 

READ(7.*)(AU  .J)  .J-10.12) 

10  CONTINUE 



REWIND  7 

CLOSE  (UNIT-7.  STATUS- ’ KEEP ' ) 

C 

800  FORMAT (5A4) 

802  FORMAT 0x.5A4) 

C 

RETURN 

END 
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SUBROUTINE  DISTR2 


*•••  PURPOSE: 

TO  DETERMINE  AIR  FLOW  DISTRIBUTION  IN  A TUBE-FINNED  COIL 
JUNE  15,1988 


••••  INPUT  DATA: 

lie 

NDEP(IIO) 

NTUBHIO.N) 

NTPS(IIO) 


1 FOR  THE  FIRST  SLAB  (-) 

2 FOR  THE  SECOND  SLAB  (-) 

NUMBER  OF  TUBE  DEPTH  ROWS  (-) 

NUMBER  OF  TUBES  IN  ROW  N IN  THE  SLAB  (-) 
NUMBER  OF  TUBES  IN  THE  SLAB  (-) 


••••  OUTPUT  DATA: 
AMR(IIO.I) 

GET(IIO.I.I) 

GET(IIO.I.2) 

IGET(IIO.I.I) 
IGET(IIO.I.2) 


RATIO  OF  AIR  MASS  FLOW  FOR  A GIVEN  TUBE  TO 
THE  TOTAL  MASS  FLOW  RATE  FOR  THE  GIVEN  SLAB  (-) 
FRACTION  OF  THE  AIR  FLOW  OF  TUBE  IGET( 1 10, 1 . 1 ) 
PASSING  THROUGH  TUBE  T. 

FRACTION  OF  THE  AIR  FLOW  OF  TUBE  IGET(I 10, I ,2) 
PASSING  THROUGH  TUBE  *1*.  IF  IGETO 10, I , 2)-0 , 
GET(II0,I.2)  IS  SET  TO  0.  (-) 

NIAIBER  OF  A TUBE  FROM  WHICH  TUBE  T IS 
GETING  AIR  (-) 

NUMBER  OF  A SECOND  TUBE  FROM  WHICH  TUBE  * I * 

IS  GETING  AIR.  IF  THE  SECOND  TUBE 
DOES  NOT  EXIST,  IGET( 1 10. 1 ,2)=0 . 


COKMON/PR I NTD/ 1 D I A . I D I A0 . NNN . 

COMyKDN/PRINT/IPR 

C0MM0N/HPHX/NSLABS.NDEP(2) .NR0W(2) .DI (2) .D0(2) .DT(2) ,TPCH(2) . 
k DPCHf2) .WIDTH(2) .FPCH(2) .FTK(2) .FMK(2) .TMK(2) ,CFM1 .BSIDE(2) . 
A NTUB(2.5)  . IFROM(2.130)  .NTPS(2)  .BSPACE^)  . 

A ACFM(2).IFIN(2).ISUR(2).SFL0W(2) 

C0MM0N/ATEST/X(2. 18)  .VX^ . 18)  .NTEST(2) 
COKWON/AIRD/IGET(2.130.2) .GET(2. 130,2) .AMR(2.130) 


•••  CHECK  AGRREMENT  BETWEEN  LOCAL  AND  TOTAL  CFM  MEASUREMENTS 
APPLY  CORRECTION  TO  THE  LOCAL  MEASUREMENT  VALUES 
CFM2-0 . 

DO  10  110=1 .NSLABS 
DO  10  I=2.NTEST(II0) 

10  CFM2=0.5*(VX(IIO.I-1)+VX(IIO.I))*(X(IIO.I)-X(IIO. 1-1 ) )*WIDTH( 1 10) 
1 +CFM2 
CFM2=CFM2*60. 

CFMR=CFM2/CFM1 
PER=ABS(CFMR-1 . )*100. 

IF(IDIA.NE.0)THEN 

WRITE(IDIA0.*)'CFM1 .CFM2=’ .CFM1 .CFM2 
WRITE(IDIA0.*) 

WRITE(IDIA0,*)’LOCAL  MEASUREMENTS  OF  AIR  DISTRIBUTION* 
IF(CFMR.GT.1  .^HEN 

WRITE(IDIA0.*) ’OVERESTIMATE  CFM  BY  ’.PER,’ 

ELSE 

WRITE(IDIA0.*) ’UNDERESTIMATE  CFM  BY  ’.PER,’  %’ 

END  IF 

WRITE(IDIA0,*) ’CORRECTION  TO  LOCAL  VELOCITY  VALUES  IS  APPLIED' 
IF(IDIA.NE.0)  WRITE(IDIA0,*) 

END  IF 

TBSIDE=BSIDE(1) 

I F ( NSLABS . EQ . 2 heS I DE=TBS I DE+BS I DE ( 2 ) 
VMEAN=CFM1/(WIDTH(1)*TBSIDE*60. ) 

IF(IDIA.NE.0)  WRITEODIA0.*)  ’VMEAN=’  .VMEAN 
DO  16  110=1 .NSLABS 
DO  16  1=1 .NTEST(IIO) 

16  VX(IIO.I)=VX(IIO.I)/CFMR 
DO  18  110=1 .NSLABS 
ACFM(IIO)-0. 

DO  18  I=2.NTEST(II0) 

18  ACFM(IIO)=0.5*(VX(IIO.I-1)+VX(IIO.I))*(X(IIO.I)-X(IIO.I-1))* 

1 WIDTH(IIO)*60.0+ACFM(IIO) 

••••  DETERMINE  AIR  DISTRIBUTION  FOR  EACH  TUBE  IN  THE  FIRST  ROW 
OF  EACH  SLAB 
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DO  70  110=1 .NSLABS 
DO  30  ITUBE=1 ,NTUB(II0,1) 

FIND  XLL  4 XRR 
IF(ITUBE.EQ. 1)THEN 
XLL=0. 

XRR=BSPACE( I IO)+0 . 5*TPCH( 1 10) 

ELSE 

XLL=BSPACE( I IO)+(REAL( ITUBE)-1 .5)*TPCH( 1 10) 
XRR=XLL+TPCH(IIO) 

IF(ITUBE.EQ.NTUB(II0.1))XRR=BSIDE(II0) 

END  IF 

00^20^ 1=1 .NTEST(IIO) 

20  IF(X(IIO, I) .GE.XR)GOTO  22 

WRITE(*.*)’BAD  INPUT  COIL  DATA.  MESSAGE  FROM  DISTR2.  LOOP  20' 
22  CONTINUE 

C****  CALC.  AVERAGE  VELOCITY  FOR  A TUBE  AND  "AMR" 

DO  25  L=1 .5 
X1=X(II0.I) 

X2-X(II0.I-1) 

VI-VX(IIO.I) 

V2=VX(II0.I-1) 

WRITE(*.225)X2.X1 .V2,V1 

225  FORMAT( ’X2,X1 .V2.V1=' .4F10.4) 

XL=AMAX1 (XLL,X2) 

WRITE(*.226)XL,XR 

226  FORMAT(*XL.XR=’ .2F8.4) 

V-V+(V1+(V2-V1)/(X2-X1)*(0.5*(XL+XR)-X1))*(XR-XL)/(XRR-XLL) 
IF(XLL.EQ.XL)GOTO  27 
XR=XL 
25  1=1-1 

WRITE(*,*) 'ERROR  IN  DISTR2.  LOOP  25.  1 10. ITUBE=' . 1 10. ITUBE 
27  CFMI=V.WIDTH(IIO)*(XRR-XLL)*60. 

AMR(IIO. ITUBE)=CFMI/ACFM(IIO) 

WRITE(*.*)’ ITUBE. V.AMR=' . ITUBE. V.AMR( 1 10. ITUBE) 

30  CONTINUE 

• •••  DETERMINE  IGET( 1 10. 1 30 .2)  4 GET( 1 10. 130 .2)  VALUES 
FOR  DEPTH  ROWS  2.3.4  AND  5 
NDEPP=NDEP(IIO) 

NTPSS=NTPS(IIO) 

C****  FIND  THE  MAX.  NUMBER  OF  TUBES  IN  ONE  DEPTH  ROW 
N1=NTUB(II0.1) 

NMAX=N1 

DO  32  ICT=2.NDEPP 
N0-NTUB(IIO.ICT) 

32  NMAX-MAX0(N0.N1  .r«t4AX) 

NROW(IIO)=NMAX 

DO  33  I-1.NTPSS 
IGET(IIO.I.2)-0 

33  GET(IIO. I .2)-0. 

N0=N1 

NL«N1 

C 

DO  50  ICT=2.NDEPP 
N=NTUB(IIO.ICT) 

NF=NL+1 

NL=NL+N 

IF(N.EQ.N0)THEN 

C****  CASE  1 . N.EO.N0 

DO  34  I-NF.NL 
IGET(IIO. I .1)=I-N0 

34  GET(IIO.I .1)=1 . 

ELSE  IF  (N.GT.N0)  THEN 

CASE  2.  N.GT.N0 

IGET(IIO.NF.1)=NF-N0 
GET(IIO.NF.1)=0. 66666 
IGET(II0.NF+1 .1)=NF-N0 
GET(II0.NF+1 .1)-0.33333 
IGET(II0.NF+1 .2)-NF-N0+1 
GET(II0.NF+1 .2)-0.5 
IGETniO.NL.1)-NF-1 
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GET(IIO.NL. 1)=0. 66666 
IGET(II0,NL-1 .1)=NF-1 
GET(II0,NL-1 .1)-0. 33333 
IGET(II0.NL-1 .2)=NF-2 
GET(II0.NL-1 .2)«=0.5 
DO  38  I=NF+2.NL-2 
IGETfIIO.I.l)=I-N0-1 
IGET(IIO.I.2)=I-N0 
GET(IIO.I.1)-0.5 
38  GET(IIO.I.2)-0.5 

ELSE 

C****  CASE  3.  N.LT.N0 

IGET(IIO.NF.1)=NF-N0 
GET(IIO.NF.1)-=1  . 

IGET(IIO.NF,2)-NF-N0+1 
GET(IIO.NF.2)-0.5 
IGET(IIO.NL.1)-NF-2 
GET(IIO.NL.l)-*0.5 
IGET(IIO.NL.2)-NF-1 
GET(IIO.NL.2)-=1 .0 
DO  40  I-NF+1.NL-1 
IGET(IIO.I.1)-=I-N0 
IGET(IIO,I.2)-=I-N0+1 
GET(IIO.I.1)-0.5 
40  GET(lIO.I.2)-0.5 

END  IF 
N0-N 

50  CONTINUE 

DO  60  I-N1+1.NTPSS 
J=IGET(IIO.I.1) 

AMR(IIO.I)»AMR(IIO.J)*GET(IIO.I .1) 

J=IGET(IIO.I.2) 

60  IF(J.NE.0)AMR(IIO. I)=AMR(IIO. I)+AMR(IIO. J)*GET(IIO.I .2) 
70  CONTINUE 
RETURN 
END 
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FUNCTION  DYNADP(P1 ,H1 . P2 , H2 , RMS . D) 

C 

C**.*  PURPOSE: 

C TO  CALCULATE  DYNAMIC  PRESSURE  DROP 

C FOR  FLOW  IN  A TUBE 

C 

C****  INPUT  DATA: 


c 

D 

- TUBE  DIAMETER 

(FT) 

c 

HI 

- REFRIG. 

ENTHALPY 

AT 

TUBE 

INLET 

c 

H2 

- REFRIG. 

ENTHALPY 

AT 

TUBE 

OUTLET 

c 

PI 

- REFRIG. 

PRESSURE 

AT 

TUBE 

INLET 

c 

P2 

- REFRIG. 

PRESSURE 

AT 

TUBE 

OUTLET 

c 

RMS 

- REFRIG. 

MASS  FLOW  RATE 

(LBM/H) 

C 

C****  OUTPUT  DATA: 

C DYNADP  - DYNAMIC  PRESSURE  DROP  (PS I) 

C 

C****  SUBPROGRAMS  CALLED  BY  DYNADP: 

C CPCV . I TRPR2 . SATPR . SATT . SATVF . VPVHS 

C 

DO100I-1 ,2 

P-P1 

H-H1 

IF(I.EQ.1)GOTO  5 

P«P2 

H-H2 

5 TG-=SATT(P) 

CALL  VPVHS(1 .TG.P.VG.HG.SG.HF) 
IF(H.LT.HF)GOTO  10 
IF(H.GT.HG)GOTO  20 
X=1 ./(HG-HF) 

X=X*(H-HF) 

VF=SATVF(TG) 

V«VF+X*(VG-VF) 

GOTO  30 

10  CPF=SATPR(5.TG) 

T-TG+(H-HF)/CPF 

V-=SATVF(T) 

GOTO  30 
20  TT-TG+10. 

CALL  CPCV(TT.P.CP.CV,GA.SO) 

T=TG+(H-HG)/CP 

PR2-H 

CALL  ITRPR2(T,2. ,P, 2. 0.001 .PR2.V.H.S) 

30  V2-V 

IF(I.EQ.1)V1-V 
100  CONTINUE 

G=0.7853982*D*D 

G=RMS/G 

G=G*G/(32 .2*144. *3600 . *3600 . ) 
DYNADP=G*(V2-V1) 

RETURN 

END 
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(BTU/LBM) 

(BTU/LBM) 

(PSIA) 

(PSIA) 
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FUNCTION  EVGUNG(RMS.X1 ,X2.TGI .HFG.VLIQ, WAP.D.AL) 

••••PURPOSE: 

TO  COMPUTE  EVAP.  HEAT  TRANSFER  COEFF.  FOR  R22  FLOW  IN  A SMOOTH 
TUBE.  JAN  13.  1989 
••••  INPUT  DATA: 

AL  - TUBE  LENGTH  (FT) 

0 - TUBE  DIAMETER  (FT) 

HFG  - REFRIG.  HEAT  OF  EVAPORATION  (BTU/LBM) 

RMS  - REFRIG.  MASS  FLOW  RATE  (LBM/H) 

TGI  - REFRIG.  SATURATION  TEMPERATURE  (F) 

VLIQ  - SPEC.  VOLUME  OF  SAT.  LIQUID  (LBM/FT^^3) 

WAP  - SPEC.  VOLUME  OF  SAT.  VAPOR  (LBM/FT^^3) 

XI  - REFRIG.  QUALITY  AT  TUBE  INLET  (-) 

X2  - REFRIG.  QUALITY  AT  TUBE  OUTLET  (-) 

••••  OUTPUT  DATA: 

EVGUNG  - EVAPORATION  GUNG  TRANSFER  COEFFICIENT  (BTU/H^F^FT*^2) 
••••  SUBPROGRAMS  CALLED  BY  EVGUNG:  SATP,  SATPR 
••••  REFERENCE: 

GUNGOR,  K.E.  AND  WINTERTON,  R.H.S..  A GENERAL  CORRELATION  FOR 
FLOW  BOILING  IN  TUBES  AND  ANNULI.  INT.J.HEAT  MASS  TRANSFER, 

VOL  29.  NO.  3.  PP.  352-358.  1986. 

DOUBLE  PRECISION  PRED.G.QFLUX.BO. FRL.XAV.PR.RELIQ.HTCLIQ.XTT. E. 

1 S.HPOOL.EVHTC 
PRED-SATP(TGI)/721 .9 
G=RMS/(0 . 7853982 *0*0) 

QFLUX=RMS*HFG*(X2-X1 )/(AL*D) 

BO=QFLUX/(G*HFG) 

FRL=(G*VLIQ)**2/(32.2*D) 

XAV=0.5*(X1+X2) 

VISLI0=SATPR(1 .TGI) 

V I SVAP=SATPR ( 2 . TG I ) 

CONLIO=SATPR(3.TGI) 

CPLIO=SATPR(5.TGI) 

PR-VISLIO*CPLIQ/CONLIQ 

RELIQ-G*(1  .-XAV)*0AISLIQ 

HTCLIO-0 . 023D0^RELIQ**0 . 8*PR**0 . 4^C0NLIQ/D 

XTT-=(1  .D0/XAV-1  .D0)**.9*(VLIQ/WAP)**.5*(VISLIQAISVAP)^^.1 

E-1  .D0+24000.D0*BO**1 .16+1 .37D0/XTT^*0.86 

I F( FRL . LT . 0 . 05D0) E-E*FRL** (0 . 1 D0-2 . D0*FRL) 

S=1  .D0/(1 .D0+1 .15D-6*E*E*RELIO**1 .17) 

I F ( FR  L . LT . 0 . 05D0 ) S-S*DSQRT (FRL) 

HPOOL-55 . D0*PRED** . 1 2* (-DLOG1 0(PRED) )••(-. 55) 

HPOOL-HPOOL*86 . 46D0* • (- . 5) • (3 . 1 525D0*QFLUX) • ^0 . 67/5 . 6745D0 

EVHTOE*HTCLIQ+S*HPOOL 

EVGUNG-EVHTC 

RETURN 

END 
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FUNCTION  EVPDP2(RMS.TGI .XI . X2 , HFG . VLIQ . WAP . VISL , VISV. AL . D) 


•••*  PURPOSE: 

TO  COMPUTE  FRICTIONAL  EVAPORATION  PRESSURE  DROP 
FOR  FLOW  IN  A TUBE  . OCT-19-87 

•••*  INPUT  DATA: 

AL  - TUBE  LENGTH  (FT) 

D - TUBE  INSIDE  DIAMETER  (FT) 

HFG  - REFRIG.  HEAT  OF  CONDENSATION  (BTU/LBM) 

RMS  - REFRIG.  MASS  FLOW  RATE  (LBM/H) 

TGI  - SAT.  TEMPERATURE  OF  REFRIG.  (F) 

XI  - REFRIG.  QUALITY  AT  TUBE  INLET  (-) 

X2  - REFRIG.  QUALITY  AT  TUBE  OUTLET  (-) 

VISL  - LIQUID  DYNAMIC  VISCOSITY  (LBM/H*FT) 

VISV  - VAPOR  DYNAMIC  VISCOSITY  (LBM/H •FT) 

VLIQ  - SPEC.  VOLUME  OF  LIQUID  REFRIG.  (FT**3/LBM) 

WAP  - SPEC.  VOLUME  OF  VAPOR  REFRIG.  (FT**3/LBM) 

••••  OUTPUT  DATA: 

EVPDP  - FRICTIONAL  EVAPORATION  PRESSURE  DROP  (PS I) 

••••  NOTE:  PIERRE  CORRELATION  USED  UNLESS  CALCULATED  PRESSURE  DROP  IS 
SMALLER  THAN  SINGLE  PHASE  PRESSURE  DROP  FOR  X*RMS  VAPOR. 

DATA  AC/1 .6654E-11/ 

G-RMS/(0 . 7853981 6*D*D) 

RE«=G*DAISL 

AKF*=778 . •HFG*  (X2-X1  )/AL 
RATIO=RE/AKF 
RATIO=AMAX1 (1 . .RATIO) 

F-0. 01 85/RAT 10* *.25 
XM=0.5*(X1+X2) 

V2PH=VLIQ+XM*(WAP-VLIQ) 

EVPDP2=AC*  F* V2PH*G*G*AL/D 
RE=G*XM*D/VISV 
F1-0.046/RE**0.2 
RATE-2.*XM**2*F1*WAP/(F*V2PH) 

IF(RATE.GT. 1 . )EVPDP2-EVPDP2*RATE 
RETURN 
END 
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SUBROUTINE  EVPHX2(KSLABS, 1 1 1 .RMASS.T1 .PI ,H1 ,ATIN, 

4 APIN.ARHIN,T2,P2.H2.X2.S2.0L) 

NOV,  22  1988 

»•••  evaporator  simulation 

•••*  THIS  PROGRAM  COMPUTES  IN  THE  FORWARD  TUBE-BY-TUBE  SCHEME 

PERFORMANCE  OF  CROSS-FLOW  AIR  HEATED  EVAPORATOR  (ONE  OR  TWO 
SLABS  ASSEMBLIES)  WITH  UP  TO  130  FINNED  TUBES  (PER  SLAB) 
**..  PLACED  IN  UP  TO  5 DEPTH  ROWS. 

••••  THE  MODEL  ACCOUNTS  FOR  AIR  AND  REFRIGERANT  DISTRIBUTION. 


INPUT  DATA: 
ACFM(IIO) 

APIN 

ARARIO(IIO.  I) 


ARHIN 

ATIN 

Dinio) 

DO(IIO) 

DPCH(IIO) 

DT(IIO) 

FLOW(IIO.M) 

FMK(IIO) 

FPCH(IIO) 

FTK(IIO) 

GET(IIO,I.1) 

GET(IIO,I .2) 


HI 

IDEPTH(IIO.M) 

IFROM(IIO,M) 

IGET(IIO,I.1) 

IGET(IIO.I.2) 

no 

III 


KSLABS 
NDEP(IIO) 
NROW(I 10) 
NTPS 

NTUBdIO.I) 

PI 

RMASS 


TPCH(IIO) 

TMK(IIO) 

T1 

WIDTH(IIO) 


AIR  MASS  FLOW  RATE  THROUGH  SLAB  (LB/H) 

AIR  INLET  PRESSURE  (PSIA) 

RATIO  OF  AIR  MASS  FLOW  RATE  FOR  A GIVEN 
TUBE  TO  THE  TOTAL  MASS  FLOW  RATE  FOR 
A GIVEN  SLAB  (-) 

AIR  INLET  RELATIVE  HUMIDITY  (-) 

AIR  INLET  TEMPERATURE  (F) 

INNER  DIAMETER  OF  TUBES  (FT) 

OUTER  DIAMETER  OF  TUBES  (FT) 

TUBE  DEPTH  PITCH  (FT) 

FIN  TIP  DIAMETER  (FT) 

FRACTION  OF  COIL  TOTAL  REFRIG.  MASS  FLOW  PASSING 
THROUGH  TUBE  M (-) 

FIN  MATERIAL  THERMAL  CONDUCTIVITY  (BTU/FT*H*F) 
FIN  PITCH  (FT) 

FIN  THICKNESS  (FT) 

FRACTION  OF  THE  AIR  FLOW  RATE  OF  TUBE 
IGET(IIO.I.I)  PASSING  THROUGH  TUBE  ’T 
FRACTION  OF  THE  AIR  FLOW  RATE  OF  TUBE 
IGET(IIO. I .2)  PASSING  THROUGH  TUBE  T 
IF  IGET(IIO.I .2)=0,  GET(II0.I.2)  IS  SET  TO  0. 
REFRIG.  ENTHALPY  AT  EVAPORATOR  INLET  (BTU/LB) 
DEPTH  ROW  OF  A TUBE  M 

NUMBER  OF  TUBE  TUBE  M RECEIVES  REFRIG.  FROM 
WHEN  COIL  WORKS  AS  EVAPORATOR  (-) 

NUMBER  OF  A TUBE  FROM  WHICH  TUBE  'I'  IS 
GETTING  AIR  (-) 

NUMBER  OF  A SECOND  TUBE  FROM  WHICH  TUBE  T 
IS  GETTING  AIR.  IF  THE  SECOND  TUBE  DOES  NOT 
EXISTS.  IGET(IIO. I .2)-0.  (-) 

SLAB  INDICATOR  FOR  THE  INPUT ED  DATA 
1 OR  2 

1 OR  2.  SPECIFIES  THE  SLAB  TO  BE  CONSIDERED 
IF  KSLABS=1.  IF  KSLABS=2.  ANY  INTEGER  VALUE 
MAY  INPUTED  FOR  1 1 1 (-) 

NUMBER  OF  EVAPORATOR  SLABS  TO  BE  SIMULATED  (-) 
NUMBER  OF  TUBE  ROW  DEPTHS  (-) 

MAX.  NUMBER  OF  TUBES  PER  DEPTH  ROW  IN  SLAB  (-) 
NUMBER  OF  TUBES  IN  THE  SLAB  (-) 

NUMBER  OF  TUBES  IN  ROW  I OF  THE  SLAB  (-) 

REFRIG.  PRESSURE  AT  EVAPORATOR  INLET  (PSIA) 
REFRIG.  MASS  FLOW  RATE  THROUGH  BOTH  SLABS 
(IF  NSLABS=2)  OR  THE  SLAB  (IF  NSLABS=1 ) 

UB/H) 

TUBE  ROW  PITCH  (FT) 

TUBE  MATERIAL  THERMAL  CONDUCTIVITY  (BTU/FT*H*F) 
REFRIG.  TEMPERATURE  AT  EVAPORATOR  INLET  (F) 

COIL  WIDTH  (FT) 


•••  OUTPUT  DATA: 
HR(IIO.I) 

H2 

PRM(II0.1 .1) 

PRM(II0,2,I) 

P2 

S2 

T2 

X2 


REFRIG.  ENTHALPY  AT  THE  OUTLET  OF  TUBE  I 
(BTU/LBM) 

REFRIG.  ENTHALPY  AT  EVAPORATOR  OUTLET  (BTU/LB) 
REFRIG.  PRESSURE  AT  I TUBE  INLET  (PSIA) 

REFRIG.  PRESSURE  AT  I TUBE  OUTLET  (PSIA) 
REFRIGERANT  PRESSURE  AT  EVAPORATOR  OUTLET  (PSIA) 
REFRIG.  ENTROPY  AT  EVAPORATOR  OUTLET  (BTU/LB*F) 
REFRIG.  TEMPERATURE  AT  EVAPORATOR  OUTLET  (F) 
REFRIG.  QUALITY  AT  EVAPORATOR  OUTLET  (-) 
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SAT.  TEMP.  OF  I TUBE  OUTLET  (FT**3/LB) 
SUBPROGRAMS  CALLED  BY  EVPHX2 : 

A I RHT3 , A I RPR . CPCV . DYNADP , F I NCON . EVGUNG , EVPDP , F I N EF2 . 1 TRPR . 

OV  LWET . SATPR . SATT . SPHDP 1 , SPHTC . UPWAO . VPVHS . VPSV . WATPR 

DOUBLE  PRECISION  DPRES 
COMMON/PR  I NTD/ 1 D I A . I D I AG0  , NNN , KfWM 
COtwWON/PRINT/IPR 
COMMON/CONST/TC , PC . VC , T FR . A J , E EP 

COMMON/HPHX/NSLABS.NDEP(2) .NROW(2) ,DI (2) .DO(2) ,DT(2) . TPCHf 2) . 

4 DPCHf2) ,WIDTH(2) ,FPCH(2) ,FTK(2) .FMK(2) ,TMK(2) .CFM1 ,BSIDE(2) . 
ft  NTUB(2.5) . IFROM(2. 130)  .NTPS(2)  .BSPACE^)  . 
ft  ACFM(2) . IFIN(2) . ISUR(2) , SFLOW(2) 

COMMON/MERG/MERGE(2,20,2) . IMER(2) , IOUT(2,20) ,NOUT(2) . 
ft  IDEPTH(2,130).FLOW(2,130).KFEED(2.130.3).KSTART(2.130).KST(2) 
COMMON/MASS/PRM (2,2,130) 

COMMON/AIRD/IGET (2, 130,2) ,GET(2, 130,2) , ARATIO(2 , 1 30) 
COMMON/OUTPR/HR (2,130) 

DIMENSION  TAIR(2,2,6),AIRN(5) ,H2IAIR(25,2),XRM(2,130) . 
ft  OMEGA(2.2,6)  ,DTHFG(2, 130) , IMC(20) , lENDO®)  .MY(130)  ,VGM(130)  , 
ftTRM(2,130) 

DIMENSION  TPIP(2,5) ,TWAT(2,5) ,HICE(2,5) ,HFGWT(2,5) ,TKICE(2,5) , 
1 KFEED0(10) 

DIMENSION  TAIRE(2, 130,2) ,WAIRE(2 , 1 30 , 2) , FSTORE(2 , 1 30) , 

1 HICEX(2,5,25) ,HFGWTX(2, 5,25) .TKICEX(2. 5,25) ,FLOW1 4(2, 130), 

2 DPM(25),WAIREX(2,130),TAEX(2,130) 

SAVE  TGI P,TAIRE,DTHFG,WAI RE. TAIR, OMEGA. MICE, HFGWT,TK ICE. TWAT, 

1 TPIP 

DATA  OMEGA/24*0./. PI/3. 1415927/ 

IF(IDIA.NE.0)  WRITE(IDIA.900)T1 ,P1 .HI .ATIN.ARHIN.RMASS 

• ESTABLISH  PARAMETERS  FOR  1 10  LOOPS 
IF(KSLABS.EQ. 1 )THEN 
112=111 
KSL=II1 
ELSE 
112=1 
KSL=2 
END  IF 

C»***  CALCULATE  INLET  PARAMETERS  FOR  REFRIGERANT  AND  AIR 
TG1=SATT(P1) 

TG1P=TG1 

CALL  VPVHS(1 .TGI .PI .VG1 .HG1 .SGI .HF1) 

X1=(H1-HF1 )/(HG1-HF1 ) 

CALL  AIRPR(1 .ATIN.APIN.ARHIN.WAIRIN.CPAIR.RAIR, 
tAMAIR.AKAIR) 

VAIR=RAIR*(TFR+ATIN)/144./APIN 
C..*.  PREPARE  OTHER  VARIABLES 
IDP=0 
HD=5000 . 

H2IAIR(1 .2)=999. 

DO  1 II0=II2.KSL 
DO  1 1=1 .NTPS(IIO) 

1 rSTORE(IIO.I)=0. 

DO  12  110=112, KSL 
IF(NNN.EQ.1.  AND  .MkW. EQ. 1 )THEN 
AAMAS=ACFM(I IO)*60./VAIR 

DTAIR=0.85*RMASS*(HG1-H1)/(CPAIR*AAMAS*NDEP(I 10)) 

DO  6 1=1 .NTPS(IIO) 

DTHFG(IIO,I)=0. 

ICT=IDEPTH(IIO.I) 

TAIRE(I10, I ,1)=ATIN-FL0AT(ICT)*DTAIR 
6 WA1RE0I0,I.1)=WAIRIN 
TAIR(II0.1 .1)=ATIN 
OMEGAfllO.I .1)=WAIRIN 
OMEGA( 110,2,1 )=WAIRIN 
DO  8 1=1 .NDEP(I 10) 

J=I  + 1 

OMEGAdlO.I  .J)=0MEGA(II0.1  .1) 

HICE(IIO. 1)=1 .E+30 
HFGWT(I 10. I)=0. 
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TKICE(IIO.  I)*=0. 

TWAT(IIO, I)=0. 

TPIP(IIO. I)=0. 

TAIR(I 10, 1 . J)=TA1R(II0.1 . I)-DTAIR 
ELSE 

DTAIR=TG1-TG1P 

DO  10  I=NTUB(II0.1)+1 .NTPS(IIO) 

10  TAIRE(IIO.I ,1)=TAIRE(I10.I ,1)+DTAIR 
END  IF 
12  CONTINUE 


••••  START  MAIN  LOOP 


DO  150  IAIR=1 ,25 

IF(IDIA.NE.0)  WRITE(IDIAG0,*) * IAIR=  ’ , lAIR 

DO  110  110=112, KSL 

ILN=0 

IF(ISUR(II0).E0.1)THEN 
HTCF1=1 .00 
HTCF2=1 .00 
DPF1=1 .00 
DPF2=1 .00 
ELSE 

HTCF1=2.00 
HTCF2=1 .45 
DPF1=0.5*HTCF1 
DPF2=0.5*HTCF2 
END  I F 

HYDD=HYDDIA(IIO) 

IIFIN=IFIN(IIO) 

NNROW=NROW(IIO) 

DDI=DI(II0I 

DDO=DO(IIO) 

DDT=DT(IIO) 

TTPCH=TPCHniO) 

DDPCH=DPCH(IIO) 

WWIDTH=WIDTH(I 10) 

FFPCH=FPCH(IIO) 

FFTK=FTK(IIO) 

FFMK=FMK(IIO) 

TTMK=TMK(I 10) 

API=PI*DDI*WWIDTH 

apo=pi*ddo*wwidth 

APM=0.5*(API+APO) 

APOAPO*  ( F FPCH-FFTK ) /F  FPCH 
AF=1 . 570796* (DDT+DDO) • (DDT-DDO) 

af=af*wwidth/ffpch 

AO=APO+AF 

AF LOW=WW I DTH • TTPCH • NNROW 

AFLOW=AFLOW* (TTPCH-DDO)/TTPCH 

AFLOW=AFLOW/FFPCH 

AFL0W=39.75*AFL0W 

WFLW=AO/TTPCH 

HP=2 . • TTMK/ ( DDO-DD I ) 

HPF=F I NOON ( DDO . FFPCH , FFTK ) 

AAMAS=ACFM(  I IO)*60.  AAIR 

AF I N=3 . 1 4 1 59 • ( DDT-DDO ) • ( DDT+DDO ) /4 . 

SEGF I N=DDT • • 3/24 . -DDT • DDT  * DDO/ 1 6 . +DDO* • 3/48 . 
SEGFIN=SEGFIN*2./(AFIN* (DDT-DDO)) 

DO  4 1=1 .NTPS(IIO) 

VGM(I)=0.0 
4 MY(I)=0 

DO  13  1=1 , IMER(IIO) 

J=MERGE(II0.I.1) 

13  MY(J)=1 
DO  15  1=1 .5 
15  AIRN(I)=0. 

C****  ASSIGN  REFRIGERANT  AND  AIR  PARAMETERS  AT  THE  INLET 
DO  110  NUM0=1 .KST(IIO) 

I=KSTART(IIO.NUMB) 

VGI=VG1 

TRI=T1 

PRI=P1 
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HRI=H1 
XRI=X1 
GOTO  32 

C.*.*  CALCULATE  PERFORMANCE  OF  NEXT  TUBE 
18  CONTINUE 

TRI=TRM0  10.  JJ) 

PRI=PRM(I 10,2. JJ) 

HR1=HR(II0,JJ) 

XRI=XRM(IIO.JJ) 

VGI=VGM(JJ) 

32  PRM(II0,1 . I)=PRI 
HRII=HRI 
XRI I=XRI 
TRII=TRI 
TRE=TRI 
PRE=PRI 
HRE=HRI 
XRE=XRI 

RMS=RMASS*FLOW( 1 10. 1 ) 

TGI=SATT(PRI) 

CALL  VPVHS(1 .TGI .PRI .VGI .HGI ,SI  .HFI) 

WAP=VG  I 
VLIO=SATVF(TGI) 

HFG=HGI-HFI 
ICT=IDEPTH(IIO. I) 

IF(ICT.EO. 1)THEN 

AMS=ARAT 1 0 ( 1 1 0 . 1 ) • AAMAS 
0MEGI=WAIRIN 
TAI=ATIN 
ELSE 

J1=IGETfII0.I,1) 

J2=IGET(IIO.I.2) 

IF(J2.EO.0)THEN 

AMS=AAMAS*ARATIO( 1 10, 1 ) 

TAI=TAIRE(II0.J1  .1) 

0MEGI=WAIRE(II0.J1  .1) 

ELSE 

AW 1 -AAMAS* AR AT  I Ofl 1 0 . J 1 ) -GET  f 1 1 0 . 1 . 1 ) 

AW2=AAMAS*ARATI0U I0.J2)*GET(II0. 1.2) 

TA1=TAIRE(II0.J1 .1) 

TA2=TAIRE(I I0.J2. 1 ) 

W1-WAIRE(II0,J1 .1) 

W2=WAIRE(II0,J2,1) 

CALL  MIXAIR(14.7,AW1 .TA1 .W1 . AW2 , TA2 .W2 . AMS . TAI .OMEGI .WMASS) 
END  IF 
END  IF 

TAE=TAIRE(II0.I.1) 

OMEGE=WAIREniO.I  .1) 

0MEGE=AMIN1 (OMEGE. OMEGI ) 

OMEAVE-0 . 5* (OMEGI+OMEGE) 

TAAV=0.5*(TAI+TAE) 

CALL  AIRPR(2.TAAV.APIN.RHA.0MEAVE.CPA.RA.AMA.AKA) 

FFFTK=FFTK+2 . -TKICEC 1 10. ICT) 

HC00=A I RHT3 ( I CT , DDO . TTPCH . DDPCH . FFPCH . FFFTK , WW I DTH . HYDD . NDEP ( 1 1 0 ) 
k IIFIN.AMS.AMA.CPA.AKA) 

HCO=HCOO* ( 1 . +HFGWT ( 1 1 0 , 1 CT ) ) 

FFEE-F I NEF2 ( TTPCH . DDPCH , DDO . FFTK . FFMK . HCO) 

UD1=HC0*(1  .-AF*0  -f^FEE^/AO) 

UD2=HICE(II0,ICT) 

ANNUL-0. 

XDRY-0 . 

CX  IF(NNN.EQ.2.0R.MMM.E0.3)THEN 

CX  WRITE(IPR,879)I .TRI .TAI .TAE.PRI ,XRI ,HRI .RMS, AMS. 

CX  1 DTHFG(IIO,JJ),Q,CPA.TAAV.HCO 

CX879  F0RMAT(1 I3,4F7.2,1F6.3.2F8.3.1F7.2.2F8.2.1F8.1 .3F8.3) 

CX  END  IF 

DP1=0. 

DP2=0 . 

AL1-0. 

IF(TAI .LT.TRI)THEN 
TSATT=SATT(PRI) 

IF(XRI . LT. 1 . )TAI-TRI+0.05 
END  IF 
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IF(I .E0.15)THEN 
WRITEf • ■ !=• . I 

WRITEf* ,•) ’PRI .TRI .XRl ,HRI .TGI .HG1=’ 

WRITE(* .‘^PRI ,TRI ,XR1 ,HRI .TGI .HGI 
END  IF 

IF(XRI . EG. 1 . )GOTO  64 
IF(XRI .GE.0.85)GOTO  54 
TUBE  SELECTION  FOR  COMPUTING  DONE 
*..•  PERFORM  HEAT  TRANSFER  4 REFRIG,  PRESSURE  DROP  CALCULATIONS 

••••  CASE  1 

••••  INLET  QUALITY  LESS  THAN  0.85 
ANNUL=1 .0 
HREE=-999 
DO  50  IHI=1 ,7 

HI=HTCF2*EVGUNG(RMS.XRI .XRE.TRI . HFG . VL IQ , WAP , DDI .WWIDTH) 

CALL  OVLWET(AO.API . APM . APO . HI . HD . HP . HPF . UD2 . UD1 . UAO . UPO . UWO) 
DPRES=UAO/(CPA*AMS) 

0=CPA.AMS*(TAI-TRI).(1 . -DEXP(-DPRES) ) 

HRE=HRI+Q/RMS 

I F(ABS(HRE-HREE) . LT . 0 . 001 )GOTO  52 
HREE=HRE 

50  XRE=(HRE-HFI)/HFG 

WRITE(IPR.*) ’50  DID  NOT  CONVERGE.  lAIR.  I .TAI  .AMS=MAIR.  1 .TAI  .AMS 
WRITE(IPR.*) ’RMS.TRI .HRI .XRI=’ .RMS.TRI .HRI ,XRI 
WRITE(IPR.*) ■HRE.XRE=’ .HRE.XRE 
52  IF(XRE. LE.0.85)GOTO  70 

••••  OUTLET  QUALITY  ABOVE  0.85  IN  A TUBE  OF  INLET  QUALITY  BELOW  0.85 
•••*  FIND  FRACTION  OF  A TUBE  WITH  QUALITY  UP  TO  0.85 
X85=0.85 

HI=HTCF2*EVGUNG(RMS.XRI ,X85.TRI .HFG.VLIQ. WAP.DDI .WWIDTH) 

CALL  OVLWET(AO.API .APM.APO.HI .HD.HP.HPF.UD2.UD1 .UAO, UPO. UWO) 

DPR  ES=UAO/ ( CPA • AMS ) 

0=CPA*AMS* (TAI-TRI ) • ( 1 .-DEXP(-DPRES) ) 

HRE=HRI40/RMS 

H85=HFI+0.85*HFG 

ANNU  L= ( H85-HR I ) • RMS/Q 

HRI=H85 

XRI=0.85 

*«•*  CASE  2 

•••*  INLET  QUALITY  AT  OR  ABOVE  0.85 

••••  CALC.  INSIDE  TUBE  H.T.C  AT  QUALITIES  0.85  AND  1.0 
54  XA=0.85 
XB=0 . 86 
HCHH=-999 . 

DO  58  IHI=1 .10 

H1A=HTCF2*EVGUNG(RMS.XA.XB.TGI .HFG.VLIQ. WAP. DDI .WWIDTH) 

CALL  OVLWET(AO.API .APM.APO.HIA.HD.HP.HPF.UD2.UD1 .UAO. UPO. UWO) 

DPR  ES=UAO/ r CPA • AMS ) 

0=CPA*AMS* (TAI-TRI ) • ( 1 .-DEXP(-DPRES) ) 

HCH=Q/RMS 

IF(ABS(HCH-HCHH).LT.0.001)GOTO  59 

HCHH=HCH 

XD=HCH/HFG 

XA=0 . 85-XD/2 . 

58  XB=0 . 85+XD/2 . 

WRITE(IPR.*) '58  DID  NOT  CONVERGE.  lAIR. I .XD=‘ . lAIR. I .XD 

59  CONTINUE 
VlSV=SATPR(2.TRn 
C0NV=SATPR(4.TRI) 

CALL  CPCV(TRI .PRI . CPR . CVR . GAR . SOR ) 

HIB=HTCF1*SPHTC(CPR.VISV. CONV . RMS .DDI) 

••••  CALC.  HEAT  TRANSFER  FOR  A TUBE  WITH  QUALITY  RANGE  0.85  - 1.00 
XRAV=XRI 
HREE=-999.0 
DO  60  IHI*=1 .12 

HI-HIA-(HIA-HIB)*(XRAV-0.85)/0. 15 

CALL  OVLWET(AO.API .APM.APO.HI .HD.HP.HPF.UD2.UD1 .UAO. UPO. UWO) 

DPR  ES*UAO/ ( CP A • AMS ) 

Q=(1 .-ANNUL)*CPA»AMS«(TAI-TRI)*(1 .-DEXP(-DPRES)) 
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HRE=HRI+Q/RMS 
XRE=(HRE-HF])/HFG 
IF(XRE.GT. 1 .0)THEN 
XRE=1 .0 

XDR Y= ( 1 . -ANNU  L ) • ( HG I -HR I ) / ( HR  E-HR  I ) 

ELSE 

XDRY=1 .-ANNUL 
END  IF 

1F(ABS(HRE-HREE) . LT. 0.001 )GOTO  63 
HREE=HRE 

60  XRAV=0.5*(XRI+XRE) 

WRITE(IPR,*) ’60  DID  NOT  CONVERGE.  lAIR.I  XDRY=  MAIR , I ,XDRY 
WRITE(IPR.*) 'XRI .XRE.HRI .HRE=’ .XRI .XRE.HRI .HRE 

63  IF(XRE. LT. 1 . )GOTO  70 
HRI=HGI 

C 

C****  CASE  3 

C****  INLET  QUALITY  EQUAL  1.;  SATURATED  OR  SUPERHEATED  VAPOR 

64  TGAV=SATT(PRI) 

VISVS=SATPR(2,TGAV) 

CONVS=SATPR ( 4 . TGAV ) 

HREE=-999. 

DO  66  IHI=1 ,7 
TRAV=0.5*(TRI+TRE) 

CALL  CPCV(TRAV,PRI . CPR , CVR , GAR . SOR) 

TT=SORT((TRAV+460. )/(TGAV+460. )) 

VISV=TT*VISVS 

CONV=TT*CONVS 

H I =HTC  F1»SPHTC(CPR.VISV. CONV . RMS . DD I ) 

CALL  OVLWET(AO,API .APM.APO.HI , HD . HP . HPF . UD2 . UD1 .UAO.UPOW.UWOW) 
DPR  ES=UAO/ ( C P A • AMS ) 

DPRES=(1 .-ANNUL-XDRY)*CPA*AMS*(1 .-DEXP(-DPRES) )/(CPR*RMS) 
0=CPR*RMS* (TAI-TRI )* ( 1 .-DEXP(-DPRES) ) 

HRE=HRI+0/RMS 
TRE=TRI+(HRE-HRI )/CPR 
IF(ABS(HREE-HRE) . LT. 0.001 )GOTO  68 
66  HREE=HRE 

WRITE(IPR.*)’52  DID  NOT  CONVERGE.  lAIR. I .HRE=’ . lAIR. I .HRE 

HEAT  TRANSFER  CALCULATIONS  FOR  TUBE  T COMPLETED 

•*••  CALCULATE  PRESSURE  DROP 
**.*  CASE  3.  QUALITY  EQUAL  1. 

68  VSP=VPSV(PRI .TRAV) 

AL1=(1 .-ANNUL-XDRY)*WWIDTH+12.*DDI 
DP1=DPF1*SPHDP1 (RMS.AL1 .DDI .VSP.VISV) 

IF(XDRY.EQ.0.)GOTO  72 
70  AL2=WWIDTH-AL1+12.*DDI 
VLIQ=SATVF(TGI) 

VISL=SATPRn  .TGI) 

VISV=SATPR(2.TGI) 

DP2=DPF2*EVPDP2(RMS.TGI .XRII .XRE.HFG.VLIQ.VGI . VISL.VISV.AL2 ,DDI ) 
72  DP3=DYNADP(PRI .HRI I .PRE.HRE.RMS.DDI ) 

PRE=PRI-DP1-DP2-DP3 

C**..  ENTHALPY  k PRESSURE  AT  TUBE  I INLET  ARE  KNOWN 
C*...  FIND  TEMP..  QUALITY  k SPEC.  VOLUME 
PRM(II0.2. I)=PRE 
TGE=SATT(PRE) 

CALL  VPVHS(1 .TGE.PCE.VGE.HGE.SGE.HFE) 

VGM( I )=VGE 
IF(HRE.LT.HGE)THEN 

XRM( 1 10. 1 )=(HRE-HFE)/(HGE-HFE) 

TRE=TGE 

ELSE 

PR0P2=HRE 

CALL  ITRPR2(TRE.2. .PRE.2.0.001 .PR0P2 . VRE .HRE .S2) 

XRM(I10. I)=1 .0 
END  IF 

TRM(IIO.I)«=TRE 
XSUP=(1 .-ANNUL-XDRY) 

HR(IIO.I)=HRE 

C 

OMECH1»0. 
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OMECH2=0 . 

C***  FIND  AIR  STATE  PAST  TUBE 
C 

C...*  REFRIGERANT  IN  ANNUAL  OR  DISPERSED  FLOW,  XSUP  < 1. 

I F(XSUP. LT. 0.9999) THEN 
I F ( XSUP . GT . 0 . 000 1 ) THEN 
0=(HGI-HRII)*RMS 
ELSE 

0=(HRE-HRII)*RMS 
END  IF 

TAE=TAI-Q/(CPA*AMS) 

TAEE=TAE+DTHFG(IIO. I) 

TTAIR=0.5*(TAI+TAEE) 

VELA=AAMAS*(460.+TTAIR)/(AFLOW*(FFPCH-FFTK-2.*TKICE(I 10. ICT))) 
TWATAI=TRI+Q/UW0 
IF(XSUP.GT. 0.00001 )THEN 
TWATAE=TWATAI 
TPIPAE=TRI+Q/UPO 
ELSE 

TWATAE=TRE+Q/UWO 
TPIPAE=(TRI+TRE)*0.5+Q/UPO 
END  IF 

CALL  AIRPRM .TWATAI .APIN.I . .OMEGWI .CPW.RW.AMW, AKW) 

CALL  AIRPRM .TWATAE.APIN. 1 . .OMEGWE.CPW.RW.AMW, AKW) 

TWATA=0 . 5* (TWATAI +TWATAE) 

OMEGW=0 . 5* (OMEGWI +OMEGWE) 

IF(OMEGI .GT.OMEGW)THEN 

DPR  ES=HCOO • APO/ (CPA  * AMS ) 

OMECHP= ( OM EG I -OM EGW ) • ( 1 . -D  EXP ( -DPR  ES ) ) 
TFM=TTAIR-FFEE*(TTAIR-TWATA) 

T END=TWAT A+ ( T FM-TWAT A ) /S  EG  F I N 
I F(TEND . GT . TAEE)TEND=TAEE 

CALL  AIRPR(1 .TEND.APIN.I . .OMEGIS.CPW.RW, AMW.AKW) 

DDTS=DDO+ ( DDT-DDO ) • ( OMEG I -OMEGW ) / ( OMEG I S-OMEGW ) 

DDTS=AMIN1 (DDTS.DDT) 

DDTS=AMAX 1 ( DDTS . DDO ) 

IF(DDTS.EQ.DDO)THEN 
OMECHF=0 . 

ELSE 

AFI N=3 . 1 41 59* (DDTS-DDO) • (DDTS+DD0)/4 . 

AFS=AF I N*  2 . •WW I DTH/FFPCH 

SEG»DDTS*  *3/24 . -DDTS*DDTS*DD0/1 6 . +DDO*  *3/48 . 
SEG=SEG«2./(AFIN. (DDTS-DDO)) 

0MEGIS=AMIN1 (OMEGIS.OMEGI ) 

OMEG FM=0MEGW+ (OMEG I S-OMEGW) *5 EG 
DPR ES=HCOO • A FS/ ( CP A • AMS ) 

0MECHF=(0MEGI-0MEGFM) • ( 1 .-DEXP(-DPRES) ) 

END  IF 

0MECH1=(1 .-XSUP)*(OMECHP+OMECHF) 

END  IF 
END  IF 

REFRIGERANT  IS  SUPERHEATED  IN  "XSUP"  FRACTION  OF  A TUBE 
IF(  XSUP. GT. 0.0001) THEN 
IF(XSUP.LT.0.9999)THEN 
0=(HRE-HGI)*RMS 
ELSE 

0=(HRE-HRII)*RMS 
END  IF 

TAE=TAI-Q/(CPA*AMS*XSUP) 

TAEE=TAE+DTHFG( 1 10. 1 ) 

IF(XSUP.LT.0.9999)THEN 
TWAT A I -TG I +Q/UWOW 

TP  I PAE«=TP  I PAE*  ( 1 . -XSUP  )+ ( ( TG I4TRE ) • . 54Q/UP0W  ) • XSUP 
ELSE 

TWAT A I =TR I +Q/UWOW 
TP I PAE= ( TR 1 1 +TRE ) « . 54Q/UP0W 
END  IF 

TWATAE=TRE+Q/UWOW 
TWATA=0 . 5* (TWATAI +TWATAE) 

CALL  AIRPRM .TWATAI .APIN.I . .OMEGWI .CPW.RW, AMW.AKW) 

CALL  AIRPRM .TWATAE.APIN. 1 . .OMEGWE.CPW.RW.AMW. AKW) 

IF(OMEGI .GT. OMEGW I) THEN 
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XCON=(OMEGI-OMEGWI )/ (OMEGWE-OMEGWI ) 

XCON=AMAX 1 f XCON . 0 . ) 

XCON=AMIN1 (XCON.I  . ) 

OMECON=AMIN1 (OMEGI .OMEGWE) 

OMEGW=0 . 5» (OMEGWI+OMECON) 

TWAT A=TWAT A I +0 . 5*  XCON • ( TWATA I -TWAT AE ) 

DPR  ES=HCOO • APO/ ( CPA • AMS ) 

OMECHP=( OMEG I -OMEGW) • ( 1 . -DEXP (-DPRES ) ) 
TFM=TTAIR-FFEE»(TTAIR-TWATA) 

T END=TWAT A+ ( T FM-TWATA) /SEGF I N 
I F(TEND . GT . TAEE)TEND=TAEE 

CALL  AIRPR(1 .TEND.APIN, 1 . .OMEGIS.CPW,RW,AMW,AKW) 
DDTS=DDO+ ( DDT-DDO ) • ( OMEG I -OMEGW  V ( OMEG I S-OMEGW ) 
DDTS=AMIN1 (DOTS, DDT) 

DDTS=AMAX1 (DDTS.DDO) 

IF(DDTS.EO.DDO)THEN 
OMECHF«=0 . 

ELSE 

AFI N=3 . 1 41 59* (DDTS-DDO) • (DDTS+DD0)/4 . 

AFS=AFIN*2. •WWIDTH/FFPCH 

SEG=DDTS*  *3/24 .-DDTS*DDTS*DDO/1 6 . +DDO**3/48 . 
SEG=SEG*2 . /(AFIN* (DDTS-DDO) ) 

OMEG I S=AM I N 1 ( OMEG I S . OMEG I ) 
OMEGFM=OMEGW+(OMEG1S-OMEGW) •SEG 
DPR  ES=HCOO*  AFS/ ( CPA • AMS ) 

OMECH  F= ( OMEG I -OMEG  FM ) * 0 ■ -D  EXP ( -DPR  ES ) ) 

END  IF 

OMECH2»=XCON*XSUP*  (OMECHP+OMECHF) 

END  IF 
END  IF 
C 

92  0MEGE-0MEGI-(0MECH140MECH2) 

TWATA=TWATAE 

tpipa=tpipae 

CALL  WATPR( TWATA. TP I PA, VELA, OMEG I .WATRO.WATK, 
tWATM . WATHFG . WATCP ) 

1F(IAIR. LT.6)THEN 

DTHFG(IIO. I)=WATHFG*(OMEGI-OMEGE)/CPA 
ELSE 

DTHFG( 1 10 . 1 )=0 . 5* (DTHFG( 1 10. 1 )+WATHFG* (OMEGI-OMEGE)/CPA) 
END  IF 

T A E=T A I - ( HR E-HR 1 1 ) • RMS/ ( CPA . AMS ) 

TAE=TAE+DTHFG(IIO. I) 

C TAE«=AM1N1  (TAE,TAI-0. 1) 

IF(AIRNf ICT) .EQ.0)GOTO  96 
AA1*=1  ./(AIRN(  ICT)+1  . ) 

AA2=AA1*AIRN(ICT) 

TAIR(II0.2. ICT+1)=AA1*TAE+AA2*TAIR(II0.2, ICT+1) 

0MEGA( 110,2, ICT+1 )=AA1 •0MEGE+AA2*0MEGA( 110,2, ICT+1 ) 
TWAT(IIO.  ICT)-AA1*TWATA+AA2*TWATniO.  ICT) 

TPIPfllO, ICT)»AA1*TPIPA+AA2*TPIP(II0, ICT) 
AIRN(ICT)=AIRN(ICT)+1 . 

GOTO  100 

96  TAIR( 1 10,2, ICT+1 )=TAE 
0MEGA( 1 10,2, ICT+1 )=OMEGE 
TWAT(IIO.ICT)=TWATA 
TPIP(II0.ICT)=TPIPA 
AIRN(ICT)=1 . 

100  TAIRE(IIO.I .2)=TAE 
WAIRE(IIO.I  .2)«=0MEGE 
C****  SELECT  A NEW  TUBE  FOR  CALCULATION 
IF(MY(I).EQ.1)THEN 
DO  102  N=2.3 
NN=KFEED(IIO.I.N) 

IF(NN.EQ.0)GOTO  104 
ILN«=ILN+1 

102  KFEED0(ILN)=NN 
END  I F 
104  JJ=I 

I-KFEED(IIO.JJ.I) 

IFM  .NE.-1)G0T0  18 
IF(ILN.NE.0)THEN 
I-KFEED0(ILN) 
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ILN=ILN-1 

JJ=IFROM(IIO.I) 

GOTO  18 
END  IF 

110  CONTINUE 

C****  ALL  TUBES  COMPUTED. 

C....  CALCULATE  REFRIGERANT  STATE  IN  THE  OUTLET  MANIFOLD  (H2  4 P2) . 
H2=0. 

P2=0 . 0 
SUMFLO0. 

DO  112  110=112, KSL 
DO  112  N=1 ,NOUT(IIO) 

I=IOUT(IIO.N) 

112  SUMFLO=SUMFLCH-FLOW(IIO.  I) 

DO  114  110=112, KSL 
DO  114  N=1 ,NOUT(IIO) 

I=IOUT(IIO,N) 

H2=FL0Wn  10,  n*HR(  1 10,  I )/SUMFL0+H2 
P2=FL0W(II0, I)*PRM(II0,2, I)/SUMFL0+P2 
CX  XOUT=XRM(I 10, I) 

CX  H0UT=HR(II0,I) 

CX  P0UT=PRM(II0,2,  I) 

CX  TOUT=TRM(IIO,I) 

CX  TSAT=SATT(POUT) 

CX  TSUP=TOUT-TSAT 

CX  WR I T E ( I PR , 864 ) I , XOUT , POUT , HOUT , TOUT , TSAT . TSUP 

CX864  F0RMAT(113,1F6.3,5F9.3) 

114  CONTINUE 

H2IAIR(IAIR,1)=H2 
DPMAX=0 . 0 

DO  118  110=112, KSL 
DO  118  N=1 ,NOUT(IIO) 

I=IOUT(IIO,N) 

DP=ABS(P2-PRM( 110,2,1)) 

IF(IDIA.NE.0)WRITE(IDIA,786)I ,IIO.DP,FL(W(IIO, I) 

786  FORMATC  1 , 1 10, DP, FLOW( 1 10, 1 )=’ ,214 ,2F7 . 3) 

118  DPMAX=AMAX1 (DPMAX.DP) 

DPM(IAIR)=DPMAX 

C****  CHECK  IF  CONVERGENCE  WAS  OBTAINED. 

IF(IAIR.EQ.1)H2S=0.0 
H2PH2— H2S‘“H2 

IF(IDIA.NE.0)  WRITE(IDIA,793)H2,H2PH2,P2,DPMAX 
793  FORMATC  H2 , H2PH2 ,P2,DPMAX  =' , 2F1 3 . 7 ,2F10 . 6) 

H2S=H2 

H2IAIR(IAIR,2)=H2PH2 

Cl 

IF(IAIR.GT.4)THEN 
IF(IAIR. LT. 15)THEN 

IF(DPM(IAIR-1).GT.0.1 .OR.DPMAX.GT.0.05)GOTO  1182 
END  IF 

IF(ABS(H2IAIR(IAIR-1 ,2)) . LT.0.3.AND.ABS(H2PH2) . LT.0. 15)GOTO  152 
END  IF 
C 



C****  CONVERGENCE  WAS  NOT  OBTAINED. 

C.**.  PREPARE  REFRIGERANT  DISTRIBUTION  FOR  THE  NEXT  LOOP. 

1182  IF(DPMAX. LT.0.05.AND. lAIR .GT . 4)THEN 
IDP=IDP+1 

DO  119  110=112, KSL 
DO  119  1=1 ,NTPS(II0) 

119  FSTORE( 1 10. 1 )=(FSTORE( 1 10. 1 )*REAL( IDP-1 )+FL0W( 1 10. 1 ) )/REAL( IDP) 
END  IF 

110=111 

IF(IAIR.LT. 15)CALL  BALFL1 (KSLABS, 1 10) 

DO  1210  110=112, KSL 
IF(IAIR.EQ.14)THEN 

DO  1992  1=1 .NTPS(IIO) 

1992  FLOW14(IIO.I)=FLOW(IIO.I) 

END  IF 

IF(IAIR.E0.15)THEN 

IF(IDIA.NE.0)WRITE(IDIA.e)'IDP=’ .IDP 
IF(IDP.GT.2)THEN 

DO  120  1=1 .NTPS(IIO) 
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120  FLOW(IIO.I)=FSTORE(IIO,I) 

ELSE 

DO  1205  1=1 ,NTPS(I 10) 

1205  FLOW(IIO.I)=0.5*(FLOW(IIO. I)+FL0W14(II0,I)) 

END  IF 
END  IF 

1210  CONTINUE 

C...*  PREPARE  AIR  SIDE  DATA  FOR  NEW  LOOP. 

C*.**  PREPARE  AIR  TEMPERATURE  AND  HUMIDITY  RATIO  FOR  EACH  TUBE. 

DO  1236  110=112, KSL 
IF(IAIR. LT. 15)THEN 

DO  121  1=1 .NTPS(IIO) 

TAIRE(IIO.I .1)=TAIREf IIO.I .2) 

121  WAIREdlO.I  .1)=WAIRE(II0.I  .2) 

IF(IAIR.EO. 12)THEN 

DO  122  1=1 .NTPS(IIO) 

TAEXdIO,  I)=0.25*TA1RE(IIO.  I .1) 

122  WAIREXdIO,I)=0.25*WAIRE(IIO.  1.1) 

END  IF 

IFdAIR.GE.13)THEN 

DO  1222  1=1 ,NTPS(I 10) 

TAEX(IIO.  I)=TAEXdIO.I)+0.25*TAIREdIO.  1.1) 

1222  WAIREX(IIO.I)=WAIREXdIO.  I)+0.25*WAIREdIO.  1,1) 

END  IF 
ELSE 

IFdAIR.E0.15)THEN 

DO  123  1=1 .NTPSCIIO) 

TAIRE(IIO.I .1)=0.25*TAIRE(IIO.I.1)+TAEX(IIO. I) 

123  WAIREdlO,  I .l)=0.25*WAIREdIO.  1 . 1 )+WAIREX(  1 10. 1 ) 

ELSE 

DO  1232  1=1 .NTPS(IIO) 

TAIREdlO,  1 .1)=0.25*TAIRE(IIO.1 .1)+0.75*TAIRE(IIO.  I .2 
1232  WAIREd  10. 1 . 1 )=0 . 25*WAIREd  10. 1 . 1 )+0 . 75*WAIREd  10. 1 .2 

END  IF 
END  IF 

1236  CONTINUE 

C***.  PREPARE  CONDESAT E/FROST  DATA;  AVERAGE  FOR  EACH  DEPTH  ROW. 

DO  149  110=112, KSL 
DO  124  1=1 .NDEP(IIO) 

J=I+1 

TAIR(II0,1 .J)=TAIR(II0.2,J) 

0MEGA(II0.1 .J)=0MEGA(II0.2.J) 

OMEGAdIO.2.J)=0. 

124  TAIR(I 10,2, J)=0. 

DO  140  1=1 .NDEP(IIO) 

TWATER=TWAT(IIO. I) 

TTAIR=TAIR(II0,1 .1) 

WWAIR=0MEGAdI0.1  .1) 

CALL  AIRPR(1 .TWATER.APIN. 1 . . WWAT ER . CCPA . RRA . 

JcAAMA.AAKA) 

TPIPE=TPIPdIO.  I) 

VELA*AAMAS* (460 .+TTAIR)/AFLOW 
VELA=VELA/(FFPCH-FFTK-2 . ‘TKICEC 1 10, 1 ) ) 

CALL  WATPR ( TWAT ER , TP I P E . VE LA . WWA I R . WATRO . WATK . 

4WATM . WATH  FG , WATCP ) 

IF(0MEGA(II0.1 .I).GT.WWATER)GOTO  125 
HICE(II0.I)=1 .E+30 
HFGWT(IIO.I)-0. 

TKICEdIO.I)=0. 

GOTO  140 

125  IF(0MEGAdI0.1.I+1).LT.0MEGAdI0.1.I))G0T0  126 
0MEGA(II0.1 .I  + 1)=0MEGA(II0.1  .1) 

HICE(IIO. I)=1 .E+30 
HFGWT(IIO.n«0. 

TKICEdIO.I)=0. 

GOTO  140 

126  WMAS=AAMAS*(0MEGA(II0.1 . 1 )-OMEGA( 1 10. 1 .1+1)) 

HFGWT (110,1 )=WATHFG* (OMEGA( 110,1,1 )-WWATER) 

CALL  AIRPR(2.TTAIR.APIN.RRRH.WWAIR,CCPA.RRA.AAMA.AAKA) 
HFGWT( 1 10. 1 )=HFGWT( 1 10, 1 )/(CCPA* (TAIR( 1 10,1,1 )-TWATER)) 
IF(TWATER.GT.32. )GOTO  128 
TK ICE( 1 10. 1 )=0 . 1 25*WMAS/(AO*NNROW*WATRO) 

TKMAX=0 . 5* ( FFPCH-FFTK ) 
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IF(TKICE(IIO. I) .GE.TKMAX)TKICE(IIO. I)=0.9*TKMAX 
GOTO  132 

128  WWW=WMAS/WFLW 

WWW=WATM • WWW/ ( WATRO • WATRO ) 

TKICE( 1 10. I ) = 1 .449E-03*WWW* *0.333 
132  HICE(IIO.I)=WATK/TKICE(IIO.I) 

IF(HICE(IIO, I) . LT.0. )HICE(1IO.I)=0. 
IF(TKICE(IIO.I).LT.0.)TKICE(IIO.I)=0. 

IF(HFGWT(IIO.I).LT.0. )HFGWT(II0. 1)=0. 

140  CONTINUE 

DO  142  ICT=1 ,NDEP(IIO) 

HICEX(IIO. ICT, IAIR)=HICE(I10. ICT) 

HFGWTXf 1 10. ICT. 1A1R)=HFGWT(II0. ICT) 

142  TKICEX(IIO. ICT.IAIR)=TKICE(IIO. ICT) 

IF(IAIR.GT.9)THEN 

DO  144  ICT=1 .NDEP(I 10) 

HICE(IIO. ICT)=0. 

HFGWTHIO.  ICT)=0. 

144  TKICE(IIO. ICT)=0. 

DO  146  ICT=1 .NDEP(I 10) 

DO  146  N=1 .4 

HICE(IIO. ICT)=HICE(IIO.ICT)+0.25*HICEX(IIO. ICT . IAIR-4+N) 
HFGWTfl 10. ICT)=HFGWTf 1 10. ICT)+0 . 25*HFGWTX( 1 10. ICT . IAIR-4+N) 
146  TKICEUIO. ICT)=TKICE(IIO. ICT)+0 . 25*TKICEX( 1 10. ICT . IAIR-4+N) 
END  IF 

149  CONTINUE 

150  CONTINUE 

* 

C****  END  OF  MAIN  LOOP 

* 

IF(IAIR.GT.25)IAIR=25 
DO  529  J=1 . lAIR 

529  IF(IDIA.NE.0)WRITE(IDIA.530)H2IAIR(J.1) .H2IAIR(J.2) 

530  FORMAT (■  H2.H2PH2=’ .2(1PE1 1 .3)) 

1-5 

DO  151  IEV-6.IAIR 

151  IF(ABS(H2IAIR(I .2)) .GT.ABS(H2IAIR(IEV.2)))I=IEV 
H2PH2=0.5*H2IAIR(I .2) 

H2=H2IAIR(I .1)+H2PH2 
IF(IDIA.NE.0)WRITE(IDIA.901)H2PH2 

152  CONTINUE 
DTT-4 . 

TOL-0.001 

CALL  PHWET2 ( P2 . H2 . DTT . TOL . T2 . V2 . S2 . X2 . TG2 ) 

TSUP2-T2-TG2 
TSUP2-AMAX 1 ( TSUP2 . 0 . ) 

QT-RMASS • ( H2-H 1 ) • SUMF  LO 
CFMIND-0. 

QL-0. 

DO  160  110=112. KSL 
CFMIND-CFMIND+ACFM( 1 10) 

160  OL-ACFM(IIO)*60.*(OMEGA(IIO.1 .1)-OMEGA(II0.1 .NDEP(II0)+1)) 

1 /(I  .-fOMEGA(II0.1  .1))  + QL 

QL=OL*WATHFG/VAIR 
OS-OT-QL 

IF(IDIA.NE.0)WRITE(1DIA.880)OT.OS.QL 
880  FORMATC  OT.OS.OL-' .3(1PE1 1 .3)) 

IF(IDIA.NE.0)WRITE(IDIA.902)T1 .PI .HI .XI .T2.P2.H2.X2.TSUP2 
C FMTON-C  FM I ND/ ( OT/ 12000.) 

IF(IDIA.NE.0)WRITE(IDIA.*) XFMIND.CFM/TON-’ .CFMIND.CFMTON 
CX  DO  310  110=112. KSL 

CX  DO  310  N=1 .NOUT(IIO) 

CX  K=IOUT(IIO.N) 

CX  PDIF-P2-PRM(II0.2.K) 

CX310  WRITE(IPR.611)P2.K.PRM(II0.2.K) .PDIF 
CX611  F0RMATfF8.3. I5.2X.F8.3.F6.3) 

900  FORMAT (/2X. ’ INPUT  DATA  TO  EVPHX2: ■//2X. ‘T’ . 

*10X. ’P’ .10X. ’H* .10X, ’TAIR’ .7X. ’RH* .9X. 'RMASS ’/eC 1PE1 1 .3)) 

901  FORMATS  EVPHX2  DOES  NOT  CONVERGE.  MAX.  ERROR-’ . F6.2. * BTU/LB’) 

902  FORMAT (/2X. ’EVAPORATOR  ITERATION: ’//2X. 

t’T’ .10X. ’P’ .10X. ’H’ .10X. ’X’ .10X. ’TSUP’ ./4(1PE11 .3)/5(1PE1 1 .3)) 
999  RETURN 

END 
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FUNCTION  FINCON(DO.FPCH,FTK) 

••••  PURPOSE:  TO  CALC.  THE  THERMAL  CONDUCTANCE  FOR  FIN-TUBE  CONTACT 
MARCH  9,  1989 

••••  INPUT:  DO  - OUTSIDE  DIAMETER  OF  TUBE  (FT) 

FPCH  - FIN  PITCH  (FT) 

FTK  - FIN  THICKNESS  (FT) 

••••  OUTPUT:  FINCON  - FIN-TUBE  CONTACT  THERMAL  CONDUCTANCE 

(BTU/(H*FTJ*F) 

••••  REFERENCE:  SHEFFIELD.  J.W..  WOOD,  R.A.  AND  SAUER.  H.J.. 

EXPERIMENTAL  I NVEST I CAT  I ON  OF  THERMAL  CONDUCTANCE 
OF  FINNED  TUBE  CONTACTS.  EXPERIMENTAL  THERMAL  AND  FLUID 
SCIENCE.  0894-1777.  1988. 

(SEE  ALSO  ASHRAE  TRANSACTION.  PAPER  NO.  3103.  1987) 

NOTE:  THIS  CORRELATION  IS  APPLICABLE  TO  MECHANICALLY  EXPANDED  COPPER 
TUBES  WITH  WITH  ALUMINIUM  FINS.  DIAMETERS  1/4  TO  5/8  INCH. 

THE  INDENTATION  DIAMATER  OF  THE  MICROHARDNESS  TEST.  D.  AND  TUBE 
EXPANSION  INTERFERENCE.  I.  ARE  GIVEN  TYPICAL  VALUES 
IN  CALCULATIONS. 

NOTE  DIMENSIONS  USED  IN  THE  CORRELATION: 

- D (MICRO  METERS) 

- ALL  OTHER  LENGTH  DIMENSIONS  IN  INCHES 

REAL  I 

DATA  D/44.2/, 1/0.0065/ 

C 

DOI=DO*12. 

FPI-1 ./(12.*FPCH) 

FTKI«=12.*FTK 

FINCON=6.902+2.889*(I*FPI*D/DOI)**0.75*(FTKI*FPI)**1 .25 
FI NCON=EXP( FINCON) 

RETURN 

END 
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FUNCT ION  F I NEF2 ( TPCH . DPCH . DR . T . AK , H ) 

••••  PURPOSE; 

CALCULATE  FIN  EFFICIENCY  USING  SCHMIDT’S  METHOD 
10-19-88 

••••  REFERENCE:  HEATING.  VENTILATING.  AND  AIR  CONDITIONING. 
MCQUISTON.F.C.  AND  PARKER.  J.D.. 

J.  WILEY  k SONS.  INC..  2-ND  EDITION.  P.478.  1982. 

••••  INPUT  DATA: 

AK  - FIN  MATERIAL  THERMAL  CONDUCTIVITY  (BTU/FT*H*F) 

DPCH  - DISTANCE  BETWEEN  TUBES  IN  NEIGHBOURING  DEPTH  ROWS 

DR  - FIN  ROOT  DIAMETER  (FT) 

H - FIN  SURFACE  HEAT  TRNSFER  COEFF.  (BTU/H*F*FT**2) 

T - FIN  THICKNESS  (FT) 

TPCH  - DISTANCE  BETWEEN  TUBES  IN  THE  SAME  DEPTH  ROW  (FT) 
••••  OUTPUT  DATA:  FINEF2  - FIN  EFFICIENCY  (-) 

REAL  L.M 

• CALC.  RATIO  ■=  REQ/R 
X-0.5*TPCH 

Y-0 . 5 * SORT ( X • X+DPCH * DPCH ) 

L-AMAXIfX.Y) 
lyNAMINI  (X.Y) 

R=0.5*DR 
PSI-M/R 
BETA=L/M 

RATIO=1 .27*PSI*SORT(BETA-0.3) 

C****  CALC.  FIN  EFFICIENCY 

FI-(RATIO-1 .)*(1 .+0.35*ALOG(RATIO)) 

M=SQRT(2.*H/(AK*T)) 

X-M*R*FI 

FINEF2-=TANH(X)/X 

RETURN 

END 
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SUBROUTINE  FRACT(NSPLIT.RN.F) 


•••*  PURPOSE: 

TO  CALCULATE  REFRIGERANT  DISTRIBUTION  AT  A SPLIT  POINT 
WITH  ONE  TUBE  SPLITING  INTO  UP  TO  20  TUBES. 

10-14-88 

INPUT; 

NSPLIT  - NUMBER  OF  TUBES  RECEIVING  REFRIG.  IN  THE  SPLIT  POINT 
RN(I)  - FLOW  RESISTANCE  (OR  # OF  TUBES)  IN  A GIVEN  CIRCUIT 
••••  OUTPUT: 

F(I)  - FRACTION  OF  REFRIGERANT  MASS  FLOW  RATE  THROUGH  A SPLIT 
FLOWING  IN  TUBE  I 

DIMENSION  RN(20).F(20).FL1(20) 

C 

IF(NSPLIT.GT.20)THEN 

WRITE(*.*)*  ERROR  IN  CALLING  FRACT,  NSPLIT  = *, NSPLIT 
RETURN 
END  IF 
C 

DO  10  1-1,20 
10  F(I)-0.0 
R1-RN(1) 

DO  20  1-2. NSPLIT 
A-R1/RN(I) 

IF(A.EQ.1.)THEN 

FL1(I)-0.5 

ELSE 

FL1(I)-1 ./(I .+A**0.571) 

END  IF 

20  FMWF(1)+(1.-FL1(I))/FL1(I) 

F(1)=1./(1.+F(1)) 

DO  30  1=2, NSPLIT 
30  F(I)-F(1)*(1 .-FL1(I))/FL1(I) 

RETURN 

END 
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FUNCTION  HYDDIA(IIO) 

• •••  PURPOSE:  TO  CALC.  HYDRAULIC  DIAMETER  FOR  SLAB  § 1 10 
JULY  15,  1988 

••••  DEFINITION: 

HYDRAQULIC  DIAMETER  -=  4*L*(CRIT.  FLOW  AREA)/(TOTAL  WETTED  AREA) 
WHERE:  L = NDEP( 1 10) •DPCH( 1 10) 

COMMON/HPHX/NS  LABS . NDEP ( 2 ) . NROW ( 2 ) . D I ( 2 ) , DO ( 2 ) . DT ( 2 ) . TPCH ( 2 ) . 
k DPCH(2).WIDTH(2).FPCH(2) .FTK(2).FMK(2) .TMK(2) .CFM1 .BSIDE(2) . 
k NTUB(2.5) . IFR0M(2. 130)  .NTPS(2)  .BSPACEU)  . 
k ACFM(2).IFIN(2).ISUR(2).SFLOW(2) 

DATA  PI/3.141592654/ 

C 

NDEPX-NDEP(IIO) 

DOX>=DO(IIO) 

TPCHX=TPCH(  no) 

DPCHX=DPCH(IIO) 

W=WIDTH(IIO) 

FTKX-FTK(IIO) 

NTPSX=NTPS(IIO) 

C CALC.  MAX.  AND  MIN.  # OF  TUBES  IN  DEPTH  ROWS 
NMAX-NTUB(IIO.I) 

NMIN-NMAX 
DO  10  I-2.5 
N*NTUB(IIO.I) 

NMAX-=MAX0(NMAX.N) 

10  NMIN=MIN0(NMIN,N) 

C CALC.  HEIGHT  OF  THE  COIL 

HGT-=REAL(NMAX-1 ) *TPCHX+0 . 75*TPCHX 
I F ( NMAX . EQ . NM I N ) HGT=HGT+0 . 5 • TPCHX 
C FLOW  AREA  BLOCKED  BY  FINS 
FN=W/FPCH(II0) 

ABFINS-FN*FTKX*HGT 

C FLOW  AREA  BLOCKED  BY  TUBES  (ADJUSTED  FOR  AREA  BLOCKED  BY  FINS) 
ABTUBS-DOX* (W-FN*FTKX) •REAL(NMAX) 

C MIN.  FLOW  AREA 

AC-W*HGT-ABF I NS-ABTUBS 

C TOTAL  HEAT  TRANSFER  AREA  ON  THE  AIR  SIDE.  ATOTAL 

AFINS=REAL(NDEPX)*DPCHX*HGT-0.25*PI*DOX*DOX*REAL(NTPSX) 

AFINS=2.*AFINS*FN 

ATUBES=PI*DOX*(W-FN*FTKX)*REAL(NTPSX) 

ATOTAL-AFINS+ATUBES 
C HYDRAULIC  DIAMETER 

HYDDIA-4*AC*NDEP( I IO)*DPCH( I IO)/ATOTAL 

RETURN 

END 
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SUBROUTINE  ITRPR2(T.DT.P. I .AER.PR0P2.V.H.S) 

••••  PURPOSE: 

TO  ITERATE  REFRIGERANT  VAPOR  PROPERTIES  FROM  GIVEN  PRESSURE 
AND  ONE  OF  THE  FOLLOWING  PROPERTIES: 

SPECIFIC  VOLUME.  ENTHALPY.  ENTROPY 

••••  INPUT  DATA: 

AER  - CONVERGENCE  PARAMETER  NOT  TO  BE  EXEEDED 
(IN  UNITS  OF  RESPECTIVE  PROPERTY) 

DT  - TEMPERATURE  INITIAL  ITERATION  STEP  (F) 

I -=1  IF  SPECIFIC  VOLUME  IS  PROP2  (-) 

-=2  IF  ENTHALPY  IS  PROP2  (-) 

* 3 IF  ENTROPY  IS  PROP2  (-) 

P - REFRIG.  VAPOR  PRESSURE  (PSIA) 

PROP2  - VALUE  OF  OTHER  THEN  PRESSURE  KNOWN  PROPERTY: 
SPECIFIC  VOLUME  (I-=1)  (FT**3/LBM) 

ENTHALPY  (1=2)  (BTU/LBM) 

ENTROPY  (1=3)  (BTU/LBM*F) 

T - APPROXIMATE  REFRIG.  VAPOR  TEMPERATURE  (F) 

••••  OUTPUT  DATA: 

H - VAPOR  ENTHALPY  (BTU/LBM) 

T - VAPOR  TEMPERATURE  (F) 

S - VAPOR  ENTROPY  (BTUABM*F) 

V - VAPOR  SPEC.  VOLUME  (FT**3/LBM) 

••••  SUBPROGRAMS  CALLED  BY  ITRPR2: 

SATT.VPVHS 

COK«wtON/PRINT/IPR 

TG-SATT(P) 

DO  100  IT-1.40 
IF(T.GT.TG)THEN 

CALL  VPVHS(2.T.P.V1 .HI  .SI  .HF) 

ELSE 

T-TG 

CALL  VPVHS(1 .T.P.V1 .HI .SI .HF) 

END  IF 

IF(I.EQ.1)PROP22=V1 
IF(I .EQ.2)PROP22-H1 
IF(I .EQ.3)PROP22=S1 
DIFF-PROP22-PROP2 
IF(ABS(DIFF).LE.AER)GOTO  110 
IF(T.EQ.TG.AND.DIFF.GT.0.)THEN 
WRITE(IPR.600)DIFF 
GOTO  110 
END  IF 

IF(IT.NE.1)GOTO  10 
T1-T 

DTT-SIGN(DT.DIFF) 

T=T-DTT 
DIFF1-DIFF 
GOTO  100 

10  DDT=(T-T1)/(DIFF-DIFF1) 

IF(ABS(DIFF1).LE.ABS(DIFF))GOTO  20 

DIFF1-DIFF 

T1-T 

20  T-T1-DDT*DIFF1 
100  CONTINUE 

WRITE(IPR.*) • ITRPR2  DID  NOT  CONVERGE.  I.P.PROP2-  ’.I.P.PROP2 
110  V-V1 
H-H1 
S-S1 
C 

600  FORMATC  ERROR  IN  CALLING  ITRPR2.  DIFF-* . 1PE1 1 .3) 

C 

RETURN 

END 
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SUBROUTINE  MIXAIR(P.AW1 .T1 .W1 ,AW2.T2.W2.AW3. r3.W3,WMASS) 


••••  PURPOSE: 

TO  CALC.  PROPERTIES  OF  THE  AIR  STREAM  RESULTED  FROM 
THE  MIXING  PROCESS  OF  TWO  WET  AIR  STREAMS 
MAR/4/1987 


«•** 


• ••* 


INPUT: 

AW1 

- 

AW2 

- 

P 

— 

T1 

- 

T2 

- 

W1 

- 

W2 

- 

OUTPUT : 

AW3 

- 

T3 

- 

W3 

- 

WMASS 

- 

MASS  FLOW  RATE  OF  THE  STREAM  #1  (LB  OF  WET 

MASS  FLOW  RATE  OF  THE  STREAM  #2  (LB  OF  WET 

AIR  TOTAL  PRESSURE  (PS I) 

TEMPERATURE  OF  THE  STREAM  #1  (F) 
TEMPERATURE  OF  THE  STREAM  #2  (F) 

HUMIDITY  RATIO  OF  THE  STREAM  #1  (LB  H20/LB 

HUMIDITY  RATIO  OF  THE  STREAM  #2  (LB  H20/LB 


AIR/H) 

AIR/H) 


DRY  AIR) 
DRY  AIR) 


MASS  FLOW  RATE  OF  THE  MIXED  STREAM  (LB  OF  WET  AIR/H) 
TEMPERATURE  OF  THE  MIXED  AIR  STREAM  (F) 

HUMIDITY  RATIO  OF  THE  MIXED  STREAM  (LB  H20/LB  DRY  AIR) 
CONDENSATION  RATE  AT  MIXING  (LB  H20/H) 


•••  NOTE:  APPLICATIONS  RANGE  OF  THIS  SUBROUTINE  IS  32<T(F)<e0. 

H(TF. W)»0.240*TF+W*( 1061 .0+0. 444*TF) 
PSATU)=EXP(0.17829*Z**3-1  .6896*2*  *2-5. 0988*2+1 3. 4353) 


WMASS<=0 . 0 
AD1-=AW1/(1  .+W1) 

AD2=AW2/(1 .+W2) 

AD3=AD1+AD2 
H1=+I(T1  ,W1) 

H2»H(T2.W2) 

HTOTAL-AD1 *H1+AD2*H2 
H3=HTOTAL/AD3 
W3*(W1 •AD1+W2*AD2)/AD3 
T3*(H3-1061 .0*W3)/(0.240+W3*0.444) 

2- 1 000 . 0/ ( T3+459 . 67 ) 

PSAT3=PSAT(2) 

WSAT3=0 . 621 98*PSAT3/(P-PSAT3) 
IF(W3.LT.(WSAT3+0.0001))  GOTO  60 


••*  MIXED  STREAM  IS  SATURATED,  CONDENSATION  DURING  MIXING  OCCURS 
T3=T3+1 .0 
DO  50  N=1 ,10 
23*1 000 . 0/(T3+459 . 67) 

PSAT3*PSAT(23) 

W3-0 . 621 98*PSAT3/(P-PSAT3) 

H3-H(T3,W3) 

WMASS-AD1 •W1+AD2*W2-AD3*W3 
HWATER— 32.01+1 .002*T3 
DH-HTOTAL-AD3*H3-WMASS*HWATER 
FR=DH/HTOTAL 

IF(ABS(DH/HTOTAL).LT. 0.00001)  GOTO  60 
IF(N.EQ.1)THEN 
DT»SIGN(0.2,DH) 

ELSE 

DT=^DH*(T31-T3)/(DH1-DH) 

END  IF 
DH1=DH 
T31=T3 
50  T3=T3+DT 

WRITE(*,*)*MIXAIR,  LOOP  50  DID  NOT  CONVERGE' 

60  AW3=AD3*(1 .+W3) 

RETURN 

END 
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SUBROUTINE  OVLWET(AO. API .APM.APO.HI .HD.HP.HPF.HL.HO, 
1 UAO.UPO.UWO) 


PURPOSE: 

TO  COMPUTE  OVERALL  HEAT  TRANSFER  COEFFICIENT 
FOR  A WET  FINNED  TUBE  (MAR  9,  1989) 

INPUT  DATA: 

AO  - TOTAL  OUTSIDE  SURFACE  AREA  (FT**2) 

API  - TUBE  INSIDE  SURFACE  AREA  (FT**2) 

APM  - SURFACE  AREA  BASED  ON  TUBE  MEAN  DIAMETER  (FT**2) 

APO  - TUBE  OUTSIDE  AREA  (FT**2) 

HD  - TUBE  INSIDE  SURFACE  DEPOSIT  HEAT  TRANSFER  COEFF. 

(BTU/H*F*FT**2) 

HI  - INSIDE  TUBE  HEAT  TRANSFER  COEFFICIENT 
(BTU/H*F*FT**2) 

HP  - TUBE  WALL  HEAT  TRANSFER  COEFFICIENT  (BTU/H*F*FT**2) 
HPF  - THERMAL  CONDUCTANCE  OF  THE  PIPE-FIN  CONTACT 
(BTU/H*F*FT**2) 

HO  - AIR-SIDE  HEAT  TRANSFER  COEFF.  FOR  WET  FINNED  TUBE 
(BTU/H*F*FT**2) 

••••  OUTPUT  DATA: 

UAO  - OVERALL  HEAT  TRANSFER  COEFF.  FOR  WET  FINNED  TUBE 
(BTU/H*F*FT**2) 

UPO  - HEAT  CONDUCTANCE  FROM  REFRIGERANT  TO  TUBE  SURFACE 
(BTU/H*F) 

UWO  - HEAT  CONDUCTANCE  FROM  REFRIGERANT  TO  WATER  (FROST) 
SURFACE  (BTU/H*F) 

U=AO/(API*HI)+AO/(API*HD)+AO/(APM*HP)+AO/(APO*HPF) 

UP(VAO/U 

U-U+1 ./HL 

UWO-AO/U 

U-U+1 ./HO 

UAO-AO/U 

RETURN 

END 
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SUBROUT I NE  PHWET2 ( P . H . DT . TOL . T . V . S . X . TG ) 

••••  PURPOSE: 

TO  FIND  REFRIGERANT  PARAMETERS 
FROM  KNOWN  PRESSURE  & ENTHALPY. 

11/5/86 

••••  INPUT: 

DT  - ITERATION  STEP  OF  TEMPERATURE  (F) 

H - REFRIG.  ENTHALPY  (BTU/LBM) 

P - REFRIG.  PRESSURE  (PSIA) 

TOL  - CONVERGENCE  TOLERANCE  OF  ENTHALPY  (BTU/LBM) 
••••  OUTPUT: 

H - REFRIG. ENTHALPY  (BTU/LBM) 

S - REFRIG.  ENTROPY  (BTU/LBM*F) 

T - REFRIG.  TEMPERATURE  (F) 

TG  - REFRIG.  SAT.  TEMPERATURE  (F) 

V - REFRIG.  SPEC  VOLUME  (FT**3/LBM) 

X - REFRIG.  QUALITY  (-) 

•••  SUBPROGRAMS  CALLED  BY  PHWET2: 

CPCV . I TRPR2 . SATT . SATVF . VPVHS 

TG-=SATT(P) 

CALL  VPVHS(1 .TG.P.VG.HG.SG.HF) 

IF(H.LT.HF)THEN 

TT-=TG 

DO  10  K=1 .5 
TAV=0.5*(TG+TT) 

T=TG- ( HF-H ) /SATPR ( 5 . TAV ) 

IF(ABS(T-TT) .LT.0.05)GOTO  11 
TT-T 

10  CONTINUE 

1 1 V=SATVF(T) 

S=SG-(HF-H)/(0 . 5* (TG+T)+459 . 67)-(HG-HF)/(459 . 67+TG) 
X-0.0 
ELSE 

IF(H.LE.HG)THEN 

T-TG 

X=(H-HF)/(HG-HF) 

V-:X*VG+0  •-X)*SATVF(TG) 

S-SG- ( HG-H ) / ( 460 . +TG ) 

ELSE 

CALL  CPCV(TG,P.CP.CV.GA.SO) 

T-TG+(H-HG)/CP 

PPOP^H 

CALL  ITRPR2(T.DT. P.2. TOL. PROP, V.H.S) 

X-1 .0001 
END  IF 
END  IF 
RETURN 
END 
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SUBROUTINE  RDATA3 


C 

C****  PURPOSE: 

C TO  READ  FROM  FILE  8 AND  PREPARE  THE  EVAPORATOR  DATA. 

C JUNE  15.  1988 

C****  NOTE: 

C AT  THE  EXIT  OF  RDATA3.  ALL  COMMON  STATEMENT  VALUES  ARE  DEFINED 

C WITH  EXCEPTION  TO  /HPHX/ACFM(2)  AND  /PRINT/IPR 

C 

COMvON/PRINT/IPR 
COMMON/RESTR/I EXP 

COMMON/HPHX/NSLABS.NDEP(2) .NROW(2) .DI (2) .DO(2) ,DT(2) .TPCH(2) , 
k DPCH(2) .WIDTH(2) .FPCH(2) .FTK(2) .FMK(2) .TMK(2) .CFMTOT . BSIDE(2) . 
ft  NTUB(2.5) . IFROM(2. 130)  .NTPS(2)  .BSPACE^)  . 
ft  ACFM(2).IFIN(2).ISUR(2).SFLOW(2) 

COMy(»ON/ATEST/X(2.18).VX(2.18)  .NTEST(2) 

DIMENSION  ATITLE(20) 

OPEN  (UNIT=8.  FILE-=*DTEV’ . STATUS=*OLD  ’ ) 

C 

C****  INPUT  EVAPORATOR  COIL  DATA 
READ(8.800)ATITLE 
WRITE(IPR.800)ATITLE 
800  FORMAT ( 20A4) 

READ(8.*)NSLABS. lEXP, CFMTOT 
C 

C****  READ  ft  REDUCE  DATA  FOR  EACH  SLAB  SEPERATELY 
DO  100  N=1.NSLABS 
READ(8.800)ATITLE 

READ(8.*)BSIDE(N) .BSPACE(N) .WIDTH(N) 

READf8.*)TPCH(N) .DPCH(N) .DI (N) .DO(N) .TMK(N) . ISUR(N) 
READ(8.*)FPCH(N) .FTK(N) .FMK(N) .IFIN(N) 

READ(8.*HNTUB(N.  I ) . 1=1  .5)  ,SFLOW(N) 

DO  12  1=1.13 
K=10*I 
9 

12  READ(8.*)(IFR0M(N.J).J^.K) 

C****  READ  AIR  DISTRIBUTION  DATA 
READ(8.*)NTEST(N) 

READ(8.*)(XfN.n.I=1.8) 

READ(8.*)(X(N, I), 1=9.16) 

READ(8.*)(VXfN.I),I=1 .8) 

READ(8.*)(VX(N.I).I-9.16) 

C****  ASSIGN  AIR  VELOCITY  AT  THE  DUCT  WALL  (IF  NOT  SPECIFIED  BY  DTEV) 
IF(X(N.1).NE.0.)THEN 
DO  14  IT=1 .NTEST(N) 

I=NTEST(N)+1-IT 

X(N.I+1)=X(N.I) 

14  VX(N.I+1)=VX(N.I) 

X(N.1)-0.0 

VX(N,1)-VX(N.2) 

NTEST(N)-NTEST(N)+1 
END  IF 
IT-NTEST(N) 

IF(X(N.IT).LT.BSIDE(N))THEN 

NTEST(N)-NTEST(N)+1 

I-NTEST(N) 

X(N.I)=BSIDE(N) 

VX(N.I)=VX(N.I-1) 

END  IF 

C****  PREPARE  DATA  FOR  CALCULATIONS 
BSIDE(N)=BSIDE(N)/12. 

BSPACE(N)=BSPACE(N)/12. 

WIDTH(N)=WIDTH(N)/12. 

TPCH(N)=TPCH(N)/12. 

DPCH(N)-DPCH(N)/12. 

DI(N)-DI(N)/12. 

DO(N)-DO(N)/12. 

FPCH(N)-FPCH(N)/12. 

FTK(N)-FTK(N)/12. 

DT(N)-S0RT(4.*TPCH(N)*DPCH(N)/3. 14159) 

DO  20  I-I.NTEST(N) 

20  X(N.I)-X(N,I)/12. 
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••••  FIND  # OF  TUBES  IN  THE  SLAB.  NTPS(N) 

FIND  # OF  DEPTH  ROWS.  NDEP(N) 

NTPS(N)-0 
DO  50  1-1 .5 

IF(NTUB(N. I) .NE.0)NDEP(N)=I 
50  NTPS(N)-NTPS(N)+NTUB(N. I) 

C****  FIND  MAX.  NUMBER  OF  TUBES  IN  A DEPTH  ROW.  NROW(N) 

M|y|Ayai0 

DO  60  ICT-1 ,NDEP(N) 

NX-NTUB(N. ICT) 

60  NMAX-MAX0(NX,NMAX) 

NROW(N)-NMAX 
100  CONTINUE 
REWIND  8 

CLOSE  (UNIT-8.  STATUS- ' KEEP ’ ) 

RETURN 

END 
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FUNCTION  SATP(TG) 


•••*  PURPOSE: 

TO  COMPUTE  REFRIGERANT  SATURATION  PRESSURE 
FOR  GIVEN  TEMPERATURE 

**•*  INPUT  DATA: 

TG  - REFRIG.  TEMPERATURE  (F) 

REFRIG.  CONSTANTS  - REFER  TO  THE  MAIN  PROGRAM. 

••••  OUTPUT  DATA: 

SATP  - REFRIG.  SATURATION  PRESSURE  (PSIA) 
COMMON/PR I NT/ I PR 

COMMON/TGPG/AG , BG . CG . DG . EG . FG . AA . BB 
COMMON/CONST/TC , PC . VC . TFR . AJ . EEP 
SAVE  TGLAST.STPLST 
DATA  TGLAST/-111 ./ 

C 

T-TG+TFR 

IF(T.LE.0.)GOTO  999 
IFfT.GT.TC)GOTO  999 
IF(ABS(TG-TGLAST) .GT.  1.0E-4)GOTO  5 
SATP-STPLST 
RETURN 

5 C-ALOG(ABSfFG-T)) 

SATP-EEP**  (AG+BGA+CG*AL0G(T)+DG*T+EG*C*  ( ( FG-T)/T)  ) 
STPLST-SATP 
TGLAST-TG 
RETURN 

999  WRITE(IPR,100) 

100  FORMAT (5X. ’ERROR  IN  CALLING  -SATP-’) 

RETURN 

END 
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FUNCTION  SATPR(I.TG) 


••••  PURPOSE: 

TO  CALCULATE  REFRIGERANT  PROPERTIES  AT  SATURATION 
JANUARY  13.  1989 

••••  INPUT  DATA: 

I -=  1 .2.3,4  OR  5 

TG  - REFRIG.  TEMPERATURE  (F) 

REFRIG.  CONSTANTS  - REFER  TO  THE  MAIN  PROGRAM. 

•••*  OUTPUT  DATA: 

SATPR(I.TG)  - SAT.  LIQUID  DYNAMIC  VISCOSITY  (LBM/H*FT) 
SATPR(2.TG)  - SAT.  VAPOR  DYNAMIC  VISCOSITY  (LBM/H*FT) 
SATPR(3.TG)  - SAT.  LIQUID  THERMAL  CONDUCTIVITY  (BTU/H*FT*F) 
SATPR(4.TG)  - SAT.  VAPOR  THERMAL  CONDUCTIVITY  (BTU/H*FT*F) 
SATPR(5.TG)  - SPEC.  HEAT  OF  SAT.  LIQUID  (BTU/LBM*F) 

DOUBLE  PRECISION  DPRES.T.AA 
COMADN/COEFPR/A (5.12) 

DIMENSION  PRLAST(5).TGLAST(5) 

DATA  TGLAST/5*-1111 ./ 

IF(ABS(TG-TGLAST(I)).GT.1 .0E-5)GOTO  5 
SATPR-PRLAST(I) 

RETURN 
5 T-TG 

DPRES-0.D0 

K— 1 

M=1 

IF(TG.GT,100.)KfHyH-6 

N-M45 

DO10J-M,N 

K-K+1 

AA“A (I  J ) 

10  DPRES»DPRES+AA*T**K 
SATPR»DPRES 
PRLAST(I)-SATPR 
TGLAST(I)-TG 
RETURN 

END 
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FUNCTION  SATT(PG) 


PURPOSE: 

TO  COMPUTE  SATURATION  TEMPERATURE  FOR  GIVEN  PRESSURE 
••••  INPUT  DATA: 

PG  - REFRIG.  PRESSURE  (PSIA) 

REFRIG.  CONSTANTS  - REFER  TO  THE  MAIN  PROGRAM. 

•••*  OUTPUT  DATA: 

SATT  - REFRIG.  SATURATION  TEMPERATURE  (F) 

COMMON/PR I NT/ I PR 

CO(«MON/TGPG/AG . BG . CG . DG . EG , FG . AA . BB 
COMMON/CONST/TC . PC , VC , TFR . AJ , EEP 
SAVE  PG LAST, SAT LST 
DATA  PG LAST/0./ 

IF(PG.LE.0.)GOTO  999 
IF(ABS(PG-PGLAST).GT.1 .0E-3)GOTO  5 
SATT-SATLST 
RETURN 

5 PLOG*ALOG(PG) 

TR-0 . 43429448*AA*PLOG+BB 

DO10ITR«=1  .30 

TRO=TR 

OALOG(ABS(FG-TRO)) 

F=AG+BG/TRO+CG*ALOG(TRO)+DG*TRO+EG* ( ( FG-TRO) •C/TRO)-PLOG 
FT— BG/TRO*  •2+CG/TROl-DG 

IF(ABS(EG) .GE.1 .E-20)FT-=FT-EG*(1 ./TR0+FG*C/TR0**2) 
TR=TRO-F/FT 

IF(ABS(TR-TRO).LE.0.05)GOTO  20 
10  CONTINUE 
999  WRITE(IPR.100)PG 

100  FORMAT (5X. ’ERROR  IN  CALLING  SATT.  PG=’ , 1PE1 1 .3) 

RETURN 

20  SATT=TR-TFR 
SATLST=SATT 
PGLAST-PG 
RETURN 
END 
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FUNCTION  SATVF(TF) 


••••  PURPOSE: 

TO  COMPUTE  SPECIFIC  VOLUME  OF  SATURATED  LIQUID  REFRIGERANT 
••••  INPUT  DATA: 

TF  - SATURATED  LIQUID  TEMPERATURE  (F) 

REFRIG.  CONSTANTS  - REFER  TO  THE  MAIN  PROGRAM. 

••••  OUTPUT  DATA: 

SATVF  - SPECIFIC  VOLUME  OF  SAT.  LIQUID  REFRIGERANT  (FT**3/LBM) 
COKMON/PRINT/IPR 

CONMON/DENSF/AL.BL.CL.DL.EL.BPL.CPL.DPL.EPL 
COMMON/CONST/TC . PC . VC . TFR . AJ . EEP 
SAVE  TF LAST. VF LAST 
DATA  TFLAST/-1111 ./ 

T-TF+TFR 

IFrT.GT.TC)GOTO  999 
IF(T.LE.0.)GOTO  999 
IF(ABS(TF-TFLAST).GT.  1.0E-3)GOTO  5 
SATVF-VLAST 
RETURN 
5 T-1.-T/TC 

SATVF=1 ./(AL+BL*T**BPL+CL*T**CPL+DL*T**DPL+EL*T**EPL) 

TFLAST-TF 

VLAST=SATVF 

RETURN 

999  WRITE(IPR.100)TF 

100  FORMAT (5X. 'ERROR  IN  CALLING  SATVF.  TF  = '.1PE10.3) 

RETURN 

END 
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FUNCTION  SPHDP(I .T.P.RMAS.AL.D) 


•••*  PURPOSE; 

TO  COMPUTE  FRICTIONAL  SINGLE  PHASE  PRESSURE  DROP  IN  A TUBE 

••••  INPUT  DATA: 

AL  - TUBE  LENGTH  (FT) 

D - TUBE  INNER  DIAMETER  (FT) 

I  •=  1 FOR  SATURATED  LIQUID  (-) 

«=  2 FOR  SATURATED  VAPOR  (-) 

«=  3 FOR  SUPERHEATED  VAPOR  (-) 

P - REFRIG.  PRESSURE  AT  INLET  (PSIA) 

RMAS  - REFRIG.  MASS  FLOW  RATE  (LBM/H) 

T - REFRIG.  TEMPERATURE  AT  INLET  (F) 

••••  OUTPUT  DATA: 

SPHDP  - PRESSURE  DROP  OVER  TUBE  LENGTH  (PSIA) 

••••  SUBPROGRAMS  CALLED  BY  SPHDP: 

SATP . SATPR . SATT . SATVF . VPSV 

IF(I.EQ.2)GOTO  1 
IF(I.EQ.3)GOTO  2 
VSP-SATVF(T) 

AMU=SATPR(1 .T) 

GOTO  3 

1 P-:SATP(T) 

2 VSP=VPSV(P.T) 

TG=SATT(P) 

AA=((T+460. )/(TG+460. ))**0.5 
AMU=AA*SATPR(2.TG) 

3 A02./(32. 174*144. *3600. **2) 

G»=0. 7853981  *0*0 

G-RMAS/G 

RE-G* D/AMU 

F-0. 046/RE* *0.2 

IF(RE.LT.2000.)F-16./RE 

SPHDP=AC*  F*VSP*AL*G*G/D 

RETURN 

END 
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FUNCTION  SPHDP1 (AM.AL.D.VSP.AMU) 

•*••  PURPOSE: 

TO  COMPUTE  FRICTIONAL  PRESSURE  DROP 
FOR  SINGLE  PHASE  FLOW  IN  A TUBE 

INPUT  DATA: 

AL  - TUBE  LENGTH  (FT) 

AM  - FLUID  MASS  FLOW  RATE  (LBM/H) 

AMU  - FLUID  DYNAMIC  VISCOSITY  (LBM/H*FT) 

D - TUBE  DIAMETER  (FT) 

VSP  - FLUID  SPECIFIC  VOLUME  (FT**3/LBM) 

ACO3.3309E-11 
G-0. 7853981 6*D*D 
G-AM/G 
RE-=G*D/AMU 
F«0. 046/RE* *0.2 
SPHDP1-ACC*F*VSP*AL*G*G/D 
RETURN 
END 
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FUNCTION  SPHTC(CP.AM.AK,RMASS.D) 

PURPOSE: 

TO  COMPUTE  SINGLE  PHASE  HEAT  TRANSFER  COEFFICIENT 
FOR  FLOW  INSIDE  TUBE 

•••*  INPUT  DATA: 

AM  - FLUID  DYNAMIC  VISCOSITY  (LBM/FT*H) 

AK  - FLUID  THERMAL  CONDUCTIVITY  (BTU/H*F*FT) 

CP  - FLUID  SPECIFIC  HEAT  AT  CONST.  PRESSURE  (BTU/LBM*F) 
D - TUBE  DIAMETER  (FT) 

RMASS  - FLUID  MASS  FLOW  RATE  (LBM/H) 

••••  OUTPUT  DATA: 

SPHTC  - SINGLE  PHASE  HEAT  TRANSFER  COEFF.  (BTU/H*F*FT**2) 

G=RMASS/ ( 0 . 7853982 • D*D ) 

RE=D*G/AM 

IF(RE.GE.2000. )GOTO10 

SPHTC»4.36*AK/D 

GOTO20 

10  PR=(AM*CP/AK)**0.4 
RE-RE**0.8 

SPHTO0 . 023*AK*PR*RE/D 
20  RETURN 
END 
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SUBROUTINE  TRACE3 


••••  THIS  PROGRAM  DETERMINES  REFRIGERANT  DISTRIBUTION  AMONG  TUBES 
IN  A CROSS  FLOW  EVAPORATOR  BASED  ON  CIRCUITRY  CONFIGURATION. 

••••  THE  EVAPORATOR  ASSEMBLY  MAY  CONSIST  OF  ONE  OR  TWO  COILS  (SLABS). 

•••*  11-18-1988 

••••  INPUT  DATA: 

NUMBER  OF  THE  TUBE  FROM  WHICH  TUBE  J RECEIVES 
REFRIGERANT.  IF  THE  TUBE  IS  CONNECTED  TO  THE 
INLET  MAINFOLD.  I FROM  IS  SET  TO  0. 

1 FOR  THE  FIRST  SLAB  (-) 

2 FOR  THE  SECOND  SLAB  (-) 

NUMBER  OF  SLABS  IN  THE  EVAPORATOR  ASSEMBLY  (-) 
NUMBER  OF  TUBES  IN  THE  ROW  N (-) 


FRACTION  OF  COIL  TOTAL  REFRIG.  MASS  FLOW 
PASSING  THROUGH  TUBE  J (-) 

DEPTH  ROW  OF  THE  TUBE  J (-) 

NUMBER  OF  SPLIT  POINTS  (-) 

NUMBER  OF  THE  TUBE  CONNECTED  TO  THE  OUTLET 
MANIFOLD.  FOUND  AS  L SUCH  TUBE  (-) 

NUMBER  OF  THE  TUBE  RECEIVING  REFRIGERANT 
FROM  TUBE  J.  FOUND  AS  N SUCH  A TUBE 
NOTE  THAT  TUBE  J CAN  FEED  UP  TO  3 TUBES 
(N  CAN  BE  1.2  AND  3).  KFEED  IS  SET  TO  -1  IF 
J TUBE  FEEDS  THE  DISCHARGE  MANIFOLD.  KFEED 
IS  SET  TO  0 IF  A TUBE  IS  NOT  FED.  (-) 

NUMBER  OF  THE  TUBE  CONNECTED  TO  THE  INLET 
MANIFOLD.  FOUND  AS  N SUCH  TUBE  (-) 

NUMBER  OF  TUBES  CONNECTED  TO  THE  INLET 
MANIFOLD  (-) 

NUMBER  OF  THE  TUBE  WHICH  FEEDS  A SPLIT  POINT. 
FOUND  AS  K SUCH  TUBE  (-) 

NUMBER  OF  TUBES  FED  BY  THE  TUBE  K (-) 

NUMBER  OF  TUBE  DEPTH  ROWS  IN  THE  SLAB  (-) 
NUMBER  OF  TUBES  CONNECTED  TO  THE  OUTLET 
MANIFOLD  (-) 

NUMBER  OF  TUBES  IN  THE  SLAB  (-) 

FRACTION  OF  TOTAL  REFRIGERANT  MASS  FLOW  RATE 
FOR  THE  COIL  FLOWING  THROUGH  SLAB  1 10 

SUBPROGRAMS  CALLED  BY  TRACE3:  FRACT 

C0MM0N/HPHX/NSLABS.NDEP(2) .NR0W(2) .DI (2) .D0(2) .DT(2) .TPCH(2) . 
k DPCH(2)  .WIDTHS)  . FPCH(2)  . FTK(2)  . FMK(2)  .TMK(2)  .CFM1  .BSIDE(2)  . 
k NTUB(2.5)  , IFR0M(2. 130)  .NTPS(2)  .BSPACEU)  . 
at  ACFM(2)  . IFIN(2)  . ISUR(2)  ,SFL0W(2) 

COMMON/MERG/MERGE(2,20.2) . IMER(2) . IOUT(2.20) .N0UT(2) . 
k IDEPTH(2. 130) .FL0W(2. 130) .KFEED(2. 130.3) .KSTART(2. 130) ,KST(2) 
DIMENSION  IDID(130) . LEFT(20) .TC(130) .RN(20) .F(20) . ITUBE(20) . 
k ISEE(20) 

C 

DO  200  II0«1.NSLABS 
C****  FIND  NUMBER  TUBES  IN  THE  SLAB 
NTPS(IIO)-=0 
DO  1 1-1.5 

IF(NTUB(IIO. I).NE.0)NDEP(IIO)-I 
1 NTPS(IIO)-NTPS(IIO)+NTUB(IIO.I) 

C 

DO  3 1-1 .NTPS(IIO) 

FLOW(IIO.I)-0. 

DO  3 J-1 .3 

3 KFEED(IIO.I.J)-0 
DO  4 1-1 .20 

4 MERGE(IIO.I.1)-0 
C 

C*»«»  FIND  TUBES  CONNECTED  TO  THE  OUTLET  MANIFOLD 
C*«**  FIND  TUBES  WHICH  FEED  SPLIT  POINTS 
IS-0 
IM-0 


IFROM(IIO.J) 


no 

NS LABS 
NTUB(IIO.N) 

OUTPUT  DATA: 
FLOW(IIO.J) 

IDEPTH(IIO.J)  - 
IMER(IIO) 
lOUT(IIO.L) 

KFEED(IIO.J.N)  - 


KSTART(IIO.N)  - 
KST(IIO) 

MERGE(IIO.K.I)  - 

MERGE(II0.K.2)  - 

NDEP(IIO) 

NOUT(IIO) 

NTPS(IIO) 

SFLOW(IIO) 
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DO  10  J=1 .NTPS(IIO) 

DO  6 1-1 .NTPS(IIO) 

IF(IFROM(IIO,I).NE.J)GOTO  6 
NM=NlvM-1 
6 CONTINUE 

IFfNM.EQ.0)GOTO  8 
IF(NM.EQ.l)GOTO  10 

TU=TU-I-1 

MERGEfllO.IM.I WJ 
MERGE(IIO.IM.2)-NM 
GOTO  10 
8 IS-IS+1 

IOUT(IIO.IS)-J 
10  CONTINUE 
NOUTf IIO)=IS 
IMER(IIO)-IM 
C 

C****  FIND  DEPTH  ROW  FOR  EACH  TUBE 
NNDEP-NDEP(IIO) 

I LAST-0 

DO  22  J-1.NNDEP 

IFIRST-ILAST+1 

I LAST- 1 LAST+NTUB ( 1 1 0 . J ) 

DO  22  I-I FIRST, I LAST 
IDEPTH(IIO.I)-J 
22  CONTINUE 
C 

C****  FIND  REFRIG.  FLOW  PATH  FROM  THE  INLET  TO  THE  OUTLET  (KFEED  ARRAY) 
DO  50  1-1 .NTPS(IIO) 

DO  48  IK-1 ,IMER(IIO) 

IF(I.NE.MERGEniO,IK.1))GOTO  48 
IDID(I)=MERGE(II0.IK.2) 

GOTO  50 
48  CONTINUE 
IDID(I)-1 
50  CONTINUE 
C 

DO  60  IS-1 .NOUT(IIO) 

I-IOUT(IIO.IS) 

KFEED(IIO.I.I)— 1 
54  J-I 

I-IFROM(IIO.J) 

IF(I.EQ.0)GOTO  60 
N-IDID(I) 

KFEED(IIO.I.N)-J 
IDID(I)-IDID(I)-1 
IF(IDID(I).EQ.0)GOTO  54 
60  CONTINUE 
C 

C****  FIND  THE  INLET  EVAPORATOR  TUBES.  KSTART( I IO,N) 

KS-0 

DO  63  1-1 ,NTPS(IIO) 
lO-IFROM(IIO.I) 

IF(IO.EQ.0)THEN 

KS-KS+1 

KSTART(IIO.KS)-I 
END  IF 
63  CONTINUE 
KST(IIO)-KS 
C 

C****  FIND  REFRIGERANT  FLOW  DISTRIBUTION 
C****  FIND  FLOW  RESISTANCE  FOR  EACH  TUBE 
DO  65  I-I.IMER(IIO) 

65  LEFT(I)-MERGE(IIO.I ,2) 

DO  70  1-1 .NTPS(IIO) 

70  TC(I)-0. 

DO  120  IL-1 ,NOUT(IIO) 

I-IOUT(IIO.IL) 

R01  . 

DO  100  IT-1 .NTPS(IIO) 

IP-IFROM(IIO,I) 

IF(IP.EQ.0)THEN 


120 


TC(I)-=RC 
GOTO  120 
END  IF 

IF(KFEED(IIO.IP.2).EQ.0)THEN 

RC-RC+1 

I-IP 

GOTO  100 
END  IF 
TC(I)-RC 

DO  75  IM-1 , IMER(IIO) 

75  IF(IP.EQ.MERGE(IIO.IM.1))GOTO  77 
77  LEFT(IM)«=LEFT(IM)-1 

IF(LEFT(IM) .GT.0)GOTO  120 
NSPLIT=MERGE( 1 10. IM.2) 

DO  90  I1»1 .NSPLIT 
N-KFEED(IIO.IP.II) 

90  RN(I1)-=TC(N) 

CALL  FRACT(NSPLIT,RN.F) 

DO  92  11-1 .NSPLIT 
N-KFEED(IIO.IP.II) 

92  FL0W(II0.N)->F(I1) 

ROTC(N)*FLOWniO.N)**1  .75 
I-IP 

100  CONTINUE 
120  CONTINUE 

NSTART-KST(IIO) 

DO  130  I-1.NSTART 
N-KSTART(IIO.I) 

130  RN(I)-TC(N) 

CALL  FRACT(NSTART.RN.F) 

DO  132  I-1.NSTART 
N-KSTART(IIO.I) 

132  FLOVIl(IIO.N)-F(I) 

C 

C****  ASSIGN  REFRIGERANT  DISTRIBUTION.  FLOW(IIO.I) 
I STORE-0 

DO  136  1-1 ,IMER(IIO) 

ITUBE(I)-0 

136  ISEE(I)-0 

DO  160  IS-1.NSTART 
I-KSTART(IIO.IS) 

IL-1 

DO  150  10-1 .NOUT(IIO) 

DO  145  IT-1 .NTPS(IIO) 

INI-KFEED(IIO.I.IL) 

IF(IN1 .EQ.-1)THEN 
IF(ISTORE.GT.0)THEN 
I-ITUBE(ISTORE) 

IL-ISEE(ISTORE) 

I STORE- I STORE-1 
GOTO  150 
END  IF 
GOTO  160 
END  IF 

IF(IL.GT.1)G0T0  137 
DO  135  11-2.3 
IN2-KFEED(II0.I.I1) 

IF(IN2.EQ.0)GOTO  137 
ISTORE-ISTORE+1 
ITUBE(ISTORE)-I 
135  ISEE(IST0RE)-I1 

137  IN2-KFEED(II0.I.2) 

IF(IN2.GT.0)THEN 

FL0W(II0.IN1)-FL0W(II0.IN1)*FL0W(II0.I) 

ELSE 

FL0W(II0.IN1)-FL0W(II0.I) 

END  IF 

I-IN1 

IL-1 

145  CONTINUE 
150  CONTINUE 
160  CONTINUE 

DO  170  1-1 .NTPS(IIO) 
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170  FLOW(IIO.I)=FLOW(IIO. I)*SFLOW(IIO) 
200  CONTINUE 
RETURN 
END 
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FUNCTION  VPSV(P.TG) 


C 

PURPOSE; 

C TO  COMPUTE  REFRIGERANT  VAPOR  SPECIFIC  VOLUME 

C****  JANUARY  11,  1989 

••••  INPUT  DATA: 

P - REFRIG.  VAPOR  PRESSURE  (PSIA) 

TG  - REFRIG.  VAPOR  TEMPERATURE  (F) 

REFRIG.  CONSTANTS  - REFER  TO  THE  MAIN  PROGRAM. 

•••*  OUTPUT  DATA; 

VPSV  - SPECIFIC  VOLUME  OF  REFRIG.  VAPOR  (FT**3/LBM) 

••••  SUBPROGRAMS  CALLED  BY  VPSV: 

SATT 

DOUBLE  PRECISION  F . FV. V. VN.Z . EMAV.T, AKTTC, ES0. ESI  . ES2 . ES3 . 
1 ES4 . ESS . ES6 , ES7 . ES32 . ES43 , ES54 . ES65 . V2 . V3 . V4 , V5 . V6 
COMMON/PR I NT/ I PR 

CO^WON/STATE/A1 .B1 .Cl .A2.B2.C2.A3,B3.C3.A4.B4.C4.A5.B5,C5, 
«tA6,B6.C6,  ALPHA.  AK 
COMMON/CONST/TC , PC . VC . TFR . AJ . EEP 
SAVE  PLAST.TGLAST.VLAST 
DATA  P LAST. TG LAST/- 1 .,-111 ./ 

C 

T-TG+TFR 

IF(T.LT.0.)GOTO999 
IF(P.LE.0.)GOTO999 
IFfABS(P-PLAST).GT.  1 .0E-4)GOTO  5 
IF(ABS(TG-TGLAST).GT.  1.0E-4)GOTO  5 
VPSV-VLAST 
RETURN 

5 TSAT-SATT(P) 

I F ( TG . LT . ( TSAT-0 . 05 ) )G0T0999 

AKTTOAK*T/TC 

ES0-=DEXP(-AKTTC) 

ES1-P 
ES2-A1 •T 

ES3-A2+B2*T+C2*ES0 

ES4«A3+B3  * T+C3  * ES0 

ES5-A4+B4*T+C4*  ES0 

ES6-A5+B5*T+C5*  ES0 

ES7-A&+B6*T4C6*ES0 

ES32-2.*ES3 

ES43-3.*ES4 

ES54-4.*ES5 

ES65-5.*ES6 

VN-=A1  •T/P 

DO  10  ITR-1.30 

V-VN 

V2-V*V 

V3-V*V2 

V4-V*V3 

V5-V*V4 

V6-V*V5 

Z-ALPHA*(V+B1) 

IF(Z.GT.150.D0)Z-150.D0 

EMAV-iDEXP(-Z) 

F-ES 1 -ES2  A-ES3  A2-ES4/V3-ES5/V4-ES6/V5-ES7  • EMAV 

FV-ES2A2+ES32A^ES43/V4+ES54/V5+ES65A6+ES7*ALPHA*EMAV 

VN-V-F/FV 

IF(DABS({VN-V)A)  • LE.  1 .D-06)GOTO20 
10  CONTINUE 
VPSV-VN+B1 
WRITE(IPR.50) 

50  FORMAT (5X. ’VPSV  DOES  NOT  CONVERGE') 

RETURN 
999  VPSVb0 

WRITE(iPR.100)TG,TSAT.P 
100  FORMAT (5X. 'ERROR  IN  CALLING  -VPSV-',  . 
ft  /'  TG.TSAT.P-' .3(1PE11 .3)) 
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RETURN 
20  VPSV=VN+B1 
VLAST-VPSV 
PLAST-P 
TGLAST-TG 
RETURN 

END 
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SUBROUTINE  VPVHS( I .TSG.P.V.H.S.HF) 

••••  PURPOSE: 

TO  COMPUTE  REFRIGERANT  PARAMETERS 
AT  AND  ABOVE  SATURATION 
*•**  JANUARY  11,  1989 

••••  INPUT  DATA: 

I - 1 FOR  SATURATED  REFRIGERANT  (-) 

-=  2 FOR  SUPERHEATED  REFRIG.  VAPOR  (-) 

P - REFRIG.  PRESSURE  (REQUIRED  FOR  1-2  ONLY)  (PSIA) 
TSG  - REFRIG.  TEMPERATURE  (F) 

REFRIG.  CONSTANTS  - REFER  TO  THE  MAIN  PROGRAM. 

•••  OUTPUT  DATA: 

H - ENTHALPY  OF  VAPOR  (BTU/LBM) 

HF  - ENTHALPY  OF  SAT.  LIQUID  (FOR  1-1  ONLY)  (BTU/LBM) 
P - SATURATION  PRESSURE  (FOR  1-1  ONLY)  (PSIA) 

S - ENTROPY  OF  VAPOR  (BTU/LBM*F) 

V - SPECIFIC  VOLUME  OF  VAPOR  (FT**3/LBM) 

•••  SUBPROGRAMS  CALLED  BY  VPVHS: 

SATP , SATT . SATVF . VPSV 

DOUBLE  PRECISION  Z.T.C.T2.T3.T4.VR.VR2.VR3,VR4.AKE.AKEXP. 

1 HI .H2.H3.H4,HFGD,HD,CD,EMAV,S1 ,S2,S3.S4,HO 
C0M40N/PRINT/IPR 

C0M40N/TGPG/AG . BG . CG . DG . EG . FG , AA . BB 

CO(yWON/STATE/A1  .B1  .Cl  .A2.B2,C2.A3.B3.C3.A4.B4.C4.A5.B5.C5. 
a:A6.B6,C6,  ALPHA,  AK 
COMMON/SPHTV/AC . BC , CC . DC , EC . FC , X , Y 
COKWON/CONST/TC . PC . VC . TFR . AJ . EEP 
SAVE  TSGLST , PLAST . VLAST , HLAST . SLAST , HFLAST 
DATA  TSGLST.PLAST/2*-111  ,/ 

C 

T-TSG+TFR 

IF(T.LE.0.D0)GOTO  999 
IF(I .EQ.2)THEN 

IF(P.LT.0,)GOTO  999 
END  IF 

IF(ABS(TSG-TSGLST).GT.  1.0E-4)GOTO  5 
IF(I .EQ.2)THEN 

IF(ABS(P-PLAST).GT.  1.0E-4)GOTO5 
END  IF 

IF(I .EQ.1)P-PLAST 
V-VLAST 
H-HLAST 
S-SLAST 
HF-HFLAST 
RETURN 

5 IF(I.EQ.1)GOTO10 
TSAT-SATT(P) 

I F ( TSG . LT . TSAT ) GOTO  999 
V«VPSV(P.TSG) 

GOTO20 

10  P-SATP(TSG) 

V-VPSV(P.TSG) 

C-D  LOG ( DABS ( FG-T )) /T 
VF-SATVF(TSG) 

HFGD-(V-VF) •P*AJ* (-BG/T+CG+OG*T--EG* ( 1 .+FG*C) ) 

20  T2-T*T 
T3-T*T2 
T4-T*T3 
VR-V-B1 

VR2-2.D0*VR*VR 
VR3-3 . D0* VR-VR2/2 . D0 
VR4-4 . D0*VR*VR3/3 . D0 
AKE-AK*T/TC 
AKEXP  -DEXP(~AKE) 

Z-ALPHA*V 

I F(Z . GT . 1 50 . D0)Z-1 50 . D0 
EMAV-DEXP(-Z) 


125 


H1-AC*T+BC*(T2/2)+CC*(T3/3.  )+DC*(T4/4.  )-F’C/T 
H2-=AJ*P*V 

H3-A2AR+A3AR2+A4/VR3+A5/VR4 
H4-C2AR+C3AR2+C4/VR3+C5ARA 
S 1 «AC  * D LOG ( T ) +BC • T+CC • T2/2 . +DC • T3/3 . -PC/  ( 2 . • T2 ) 

52- =AJ*A1*DLOG(VR) 

53- B2/VR+B3/VR2+B4/VR3+B5/VR4 

54- H4 

IF(ABS(A6).LE.1 .E-20)GOTO30 
H0-EMAV 

I F ( ABS ( C 1 ) . GT . 1 . E-20 ) H0-H0-C 1 • D LOG ( 1 . D0+EMAV/C 1 ) 

H0-H0/ALPHA 

H3-H3+A6*H0 

H4»H4-C6*H0 

S3=S3+B6*H0 

S4=S4-C6*H0 

30  HD-H1+H2+AJ*H3+AJ*AKEXP*(1 .+AKE)*H4+X 
S-=S  1 +S2-A  J * S3+A  J • AK  EXP  • AK/TC  * S4+Y 
HFG-HFGD 
H-HD 

IF(I .EQ.1)HF-H-HFG 
TSGLST-TSG 
PLAST-P 
VLAST-V 
HLAST-H 
SLAST-S 
HFLAST-HF 
RETURN 

999  WRITE(IPR.100)TSG.P 

100  FORMAT (5X, ’ERROR  IN  CALLING  -VPVHS-’ ,2(1PE11 .3)) 
RETURN 

END 
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SUBROUT I NE  WATPR ( TW . TP , VA . WA . WATRO . WATK . WATM . WATHFG . WATCP ) 
••••PURPOSE : 

TO  COMPUTE  WATER  AND  FROST  PROPERTIES 
••••  INPUT  DATA: 

TW  - WATER  (FROST)  TEMPERATURE  (F) 

TP  - TUBE  TEMPERATURE  (F) 

WA  - AIR  HUMIDITY  RATIO  UbM  H20/LBM  DRY  AIR) 

VA  - AIR  VELOCITY  (FT/SEC) 

••••  OUTPUT  DATA: 

WATRO  - DENSITY  OF  WATER  (FROST)  (LBM/FT^^3) 

WATK  - THERMAL  CONDUCTIVITY  OF  WATER  (FROST)  (BTU/H^F^FT) 
WATM  - DYNAMIC  VISCOSITY  OF  WATER  (FROST)  (LBM/H^FT) 
WATHFG  - WATER  HEAT  OF  CONDENSATION  OR 

FROST  HEAT  OF  SUBLIMATION  (BTU/LBM) 

WATCP  - SPECIFIC  HEAT  OF  WATER  (FROST)  (BTU/LBM^F) 

DIMENSION  ARO(5) .AK(5) .AM(5) .AHFG(5) 

ARO(1)-0.11647E03 

AROf2)— 0.40054E00 

ARO(3)-0.10815E-02 

ARO(4)->-0.12387E-05 

ARO(5)«0.49002E-09 

AK(1)— 0.27694E00 

AK(2)-0.45215E-03 

AK{3)=0.49008E-05 

AK(4)-^0.88613E-08 

AK(5)-=0.41387E-11 

AMM)=0.79424E03 

AM(2)— 0.47589E01 

AM(3)-0.10622E-01 

AM(4)— 0.10416E-04 

AM(5)-0.37690E-08 

AHFG(1)-0.31514E04 

AHFG(2)«»-0.13714E02 

AHFG(3)=0.35945E-01 

AHFG(4)->-0 . 43525E-04 

AHFG(5)-0. 19695E-07 

TWR=TW+460 . 

TPR*TP+460. 

WATRO-0 . 

WATK=0 . 

WATk^0 . 

WATHFG=0 . 

IF(TW.LE.32.)G0T0  100 

DO10I=1 ,5 

J=I-1 

WATRO=WATRO+ARO( I ) •TWR^ • J 
WATK=WATK+AK ( 1) •TWR • • J 
WATI^WATM+AM(  I ^•TWR^^J 
10  WATHFG=WATHFG+AHFG(I)*TWR^^J 
WATCP=1 . 

GOTO  200 

100  B1=-11 .9521+0.02422^TPR+35.5498^WA 

ft-9 . 1742 E-07.VA+3. 1 138E-09^VA^TPR-0. 03838 
B2=(  13. 1606-0. 021 33^TPR-81 .955^WA)/(32.018-TP) 

WATRO=10  *EXP(B1+B2^(TW-TP)) 

WATK=0 . 01 2138+3.8909E-03^WATROf5. 1409E-06^WATRO^^3 

WATM=1 . E25 
WATHFG=1219 
WATCP=0.46 
200  RETURN 
END 
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APPENDIX  C.  Example  of  Run  of  the  Program,  EVSIM 


Below  is  a printout  from  execution  of  EVSIM.  The  simulated  coil  is  shown 
in  Figure  Al . The  coil  data  file  is  presented  in  Table  A4 . The  working 
fluid  is  Refrigerant  22. 


***************************************************************************** 
EVSIM  (VER.  1.1):  SIMULATION  OF  AN  EVAPORATOR  COIL 

***********************************************')t***'jSr'}t'*')!r*'3!r')!r*'!!r'*'!Sr 

March  14,  1989 
COIL  INFORMATION: 

***DTEV***  A- SHAPE  COIL,  3 DEPTH  ROWS,  16  TUBES  PER  ROW. 


NUMBERS  OF  SLABS  IN  THE  ASSEMBLY:  2 

NUMBER  OF  EXPANSION  DEVICES:  2 

TOTAL  CFM:  1120 

AIR  CONDITION: 

AIR  DRY  BULB  TEMPERATURE:  80.00 

AIR  RELATIVE  HUMIDITY:  0.51 

REFRIGERANT  22 
REFRIGERANT  CONDITIONS: 

REFRIG.  QUALITY  AT  INLET:  0.20 

REFRIG.  SAT.  TEMPERATURE  AT  EXIT:  45.00 
REFRIG.  SUPERHEAT  AT  EXIT:  10.00 


F 


F 

F 
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SIMULATION  RESULTS 

*RESULTS  FROM  THE  LAST  ITERATION  LOOP: 

TOTAL  CAPACITY:  34426.  BTU/H 
LATENT  CAPACITY:  9273.  BTU/H 
REFRIGERANT  MASS  FLOW  RATE:  494.2  LB/H 
REFRIGERANT  PARAMETERS  AT  THE  COIL  OUTLET: 
SATURATION  TEMPERATURE:  45.0  F 
SUPERHEAT:  10.1  F 
PRESSURE:  90.7  PSIA 
ENTHALPY:  110.4  BTU/LB 

INFORMATION  ON  REFRIGERANT  LEAVING  OUTLET  TUBES 


SLAB  # 

TUBE  # 

P 

T 

TSUP 

X 

RMS 

(-) 

(-) 

(PSIA) 

(F) 

(F) 

(-) 

(LB/H) 

1 

1 

90.70 

66.6 

21.6 

1.000 

123.46 

1 

16 

90.70 

63.8 

18.8 

1.000 

123.66 

2 

1 

90.68 

64.6 

19.6 

1.000 

120.64 

2 

16 

90.68 

45.0 

0.0 

0.962 

126.48 
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★RESULTS  FOR  THE  REQUESTED  REFRIGERANT  SUPERHEAT:  10.0 
(INTERPOLATED  FROM  RESULTS  OF  LAST  TWO  ITERATIONS) 


TOTAL  CAPACITY:  34450.  BTU/H 
LATENT  CAPACITY:  9291.  BTU/H 
REFRIGERANT  MASS  FLOW  RATE:  494.8  LB/H 
REFRIGERANT  PARAMETERS  AT  THE  COIL  OUTLET: 


SATURATION  TEMPERATURE 
SUPERHEAT 
PRESSURE 
ENTHALPY 


45.0  F 

10.0  F 
90.7  PSIA 

110.3  BTU/LB 
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